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1.  Problem  Statement 


The  aim  of  this  research  work  is  to  build  and  test  a  large  area ,  high  resolution ,  digital  imaging  system, 
which  can  be  used  in  the  Army  and  in  medical  care.  Digital  images  are  obtained  from  CCD  sensor  arrays. 
CCD  imaging  has  many  advantages  over  traditional  film  photography,  such  as  its  linearity,  larger  dynamic 
range  and  higher  sensitivity.  Despite  all  of  these  advantages,  however,  one  of  its  major  shortcomings  is  the 
low  spatial  resolution.  “Butting”  technology  is  an  economic  solution  for  this  problem,  where  one  or  more 
CCD  subarrays  are  butted  together  to  achieve  large  area  imaging  (Fig.  1). 

missing  information  about  50  micron  width 


Figure  1 .  Demonstration  of  two  scintillator/fiber/CCDs  butted  together. 

Problem  with  butting  is  that  no  matter  how  close  the  two  subarrays  are,  there  will  always  be  a  gap  on  the 
order  of  one  to  two  pixels  which  is  about  50  micron  wide.  Therefore,  optimal  estimation  algorithms  are 
required  to  restore  the  missing  data  while  at  the  same  time  recover  the  image  from  other  defections,  such  as 
blur,  noise,  radiometric  and  geometric  distortions. 


Figure  2.  An  original  captured  image  from  our  imaging  system. 

Fig.  2  shows  an  original  image  captured  from  our  imaging  system.  Severe  vignetting-type  radiometric  dis¬ 
tortion  and  pincushion-type  geometric  distortion  can  be  observed.  The  problems  we  need  to  solve  include: 

•  imaging  system  assembly; 

•  large  area  image  display; 

•  radiometric  correction; 

•  geometric  correction;  and 

•  missing  data  estimation  along  with  deblurring  and  denoising. 
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2.  Results  Summary 


Five  modules  are  designed  and  implemented  in  this  system:  image  acquisition,  image  display,  image 
enhancement,  image  correction,  and  image  restoration  (Fig.  3). 

/ - " 

Image  Acquisition 
Module 

v _ ~  _ , 

5 


Figure  3.  System  module  achievements. 

The  image  acquisition  module  involves  a  1  x  2  CCD/scintillator/fiber  optic  taper  combinations,  the  CCD 
driver  assembly,  the  A/D  converter,  and  the  control  circuits. 

In  our  system,  the  spatial  resolution  of  a  complete  image  is  roughly  1152  x  2484,  involving  alx21152x 
1242  arrays.  Since  most  of  the  current  commercial  monitors  can  only  display  image  with  dimensions  1024  x 
1024  or  smaller,  how  to  efficiently  and  conveniently  display  larger  images  on  smaller  display  planes  is  a 
problem  which  needs  to  be  solved. 

The  image  correction  module  consists  of  radiometric  and  geometric  corrections.  Both  radiometric  and  geo¬ 
metric  distortions  are  caused  by  the  distorted  shape  of  the  fiber  optic  tapers.  We  use  polynomial  approxima¬ 
tion  to  model  a  series  of  concentric  ellipses  formed  by  vignetting-type  radiometric  distortion.  We  use  both 
polynomial  approximation  and  thin-plate  spline  interpolation  to  correct  the  geometric  distortion  (mainly 
pincushion  distortion).  We  show  that  thin-plate  spline  interpolation  performs  better  than  polynomial  interpo¬ 
lation  since  it  maps  all  the  control  points  to  their  correspondence  exactly. 

The  image  restoration  module  is  used  to  estimate  missing  data  along  with  deblurring  and  denoising.  We  use 
two  approaches  to  estimate  missing  data  in  a  blurred,  noisy  image:  the  “consistency  method  with  separable 
deblurring”,  and  the  “maximum  a-posteriori  probability  (MAP)  estimate  with  mean  field  annealing  (MFA)”. 
Both  methods  try  to  convert  an  ill-posed  image  restoration  problem  to  a  well-posed  one  by  relaxation  theory. 
The  consistency  method  is  able  to  recover  missing  data  exactly  from  a  blurred  image  based  on  a  few  assump¬ 
tions.  However,  this  exact  restoration  ability  is  limited  by  these  assumptions,  and  the  solution  is  unstable  for 
large  noise  and  inaccurate  estimation  of  the  point  spread  function  (PSF).  The  MFA  approach  solves  the 
MAP  estimate  by  global  optimization  technique.  It  restores  the  missing  data  not  so  exact  as  the  consistency 
method,  but  it  is  more  stable  to  large  noise  and  inaccurate  estimation  of  the  PSF. 

Experiments  are  carried  out  on  both  synthetic  data  and  real  data.  The  experimental  results  show  that  the 
imaging  system  we  developed  exhibits  good  performance,  and  the  algorithms  high  efficiency. 

2.1  Image  Acquisition  and  Display 

Image  acquisition  is  the  first  step  in  any  imaging  systems.  For  CCD  imagers,  the  image  acquisition  module 
functions  as  a  combination  of  photo-detector,  amplifier,  sampler,  and  quantizer.  The  image  display  module 


Image  Correction  Module 
(Radiometric  and  Geometric  Correction) 
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reformulates  the  stream  of  sampled  and  quantized  digital  data  to  a  2-D/3-D  matrix,  and  paints  it  on  a  2-D 
image  plane/computer  monitor.  The  whole  process  is  illustrated  in  Fig.  4. 


Source 


With  a  1MHz  readout  rate,  1M  sample  rate,  1152  vertical  pixels  and  1242  horizontal  pixels  (equating  to 
1,430,784  total  pixels  in  an  array),  the  performance  of  our  data  acquisition  module  can  be  summarized  as 
follows:  1.43s  minimum  is  required  to  readout  a  complete  frame;  1.24ms  minimum  is  required  to  readout  a 
complete  line;  and  l|is  is  required  to  readout  1  pixel  of  information.  The  detail  design  of  image  acquisition 
module  is  described  in  Appendix  7. 1 . 

In  image  display  module,  we  use  a  secondary  display  board  combined  with  the  host  video  card  to  achieve  the 
so-called  “dual-screen  mode  loupe”  display.  This  technique  can  display  complete  images  on  host  monitor, 
and  segments  of  images  on  the  secondary  monitor  for  in-depth  analysis,  such  as  zooming,  autoscaling,  con¬ 
trast  stretching,  pseudo-color  enhancement,  and  missing  data  estimation.  The  design  idea  is  explained  in 
Appendix  7.2. 

2.2  Image  Correction 

Based  on  the  analysis  of  sources  of  radiometric  and  geometric  distortions  (Appendix  7.3),  we  implemented 
polynomial  approximation  and  thin-plate  spline  interpolation  to  do  the  correction. 

2.2.1  Radiometric  Correction 

Radiometric  distortion  is  modeled  by  a  series  of  concentric  ellipses.  The  conventional  polynomial  approxi¬ 
mation  is  modified  to  be  able  to  inverse  this  model  (See  details  in  Appendix  7.4).  Fig.  5  shows  two  real 
images  (a  grid  template  and  a  fiat  frame)  captured  from  our  imaging  system,  and  the  corrected  results.  It  is 
apparent  that  the  concentric  ellipse  effect  disappeared  in  the  corrected  images. 

Table  1  is  a  rough  measurement  on  how  the  correction  algorithm  works.  Since  the  correction  algorithm  is  to 
correct  the  degraded  brightness  with  central  pixel  brightness  as  the  reference,  the  mean  brightness  of  the  cor¬ 
rected  image  should  be  always  larger  than  that  of  the  measured  image.  On  the  other  hand,  for  the  grid  tem¬ 
plate,  the  corrected  image  show  better  contrast  between  white  pixel  areas  and  dark  grid  lines,  therefore,  the 
variance  of  the  corrected  image  is  larger  than  that  of  the  measured  image;  while  for  the  flat  frame  image,  the 
variance  is  much  less  than  that  of  the  measured  image  since  the  corrected  image  has  a  homogeneous  bright¬ 
ness  distribution. 
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TABLE  1.  Mean  and  Variance  comparison. 


■HUSH 

grid  template 

measured  image 

147.24 

hhIhsh 

corrected  image 

183.57 

9601.04 

flat  frame 

measured  image 

188.02 

285.42 

corrected  image 

251.60 

52.67 

Figure  5.  Results  from  radiometric  correction  of  a  grid  template  (left)  and  a  flat  frame  (right). 


We  also  use  histograms  and  profiles  to  evaluate  the  correction  results.  Fig.  6  shows  histograms  of  measured 
images  and  corrected  ones.  We  can  see  that  the  histogram  of  the  measured  grid  template  has  two  peaks:  one 
is  from  dark  grid  lines;  the  other  smoother  peak  is  from  white  pixel  areas  which  shows  the  shape  of  a  Gaus¬ 
sian  because  of  radiometric  distortion.  The  histogram  of  the  corrected  grid  template  also  has  two  peaks:  one 
is  from  dark  grid  lines;  the  other  peak  is  an  impulse  (Dirac  delta  function)  instead  of  a  Gaussian,  which  indi¬ 
cates  the  successful  correction.  Histograms  of  the  measured  flat  frame  and  the  corrected  frame  show  similar 
characteristics  as  those  of  the  grid  template  except  that  the  smaller  peak  is  from  the  boundary  of  the  image 
area. 
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grid  template  flat  frame 


Figure  6.  Comparison  of  histograms  of  the  measured  images  (top),  and  the  corrected  images  (bottom).  In  each  figure, 
the  left  plot  is  for  the  image  from  the  left  sensor,  and  the  right  plot  is  that  from  the  right  sensor. 

Fig.  7  is  a  set  of  profiles  used  to  measure  performance  of  the  algorithm  in  further  detail.  These  profiles  are 
from  the  same  columns/rows  of  the  measured  grid  template  image  and  the  corrected  image.  Profiles  on  the 
left  column  are  from  the  measured  image  which  show  a  Gaussian  envelope;  while  the  profiles  to  the  right  are 
from  the  corrected  image  that  show  a  line  envelope. 
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Profiles  of  the  measured  image.  Profiles  of  the  corrected  image. 


Figure  7.  Comparison  of  profiles  from  the  same  columns  and  row  of  the  measured  grid  template  image  and  the 
corrected  image. 


2.2.2  Geometric  Correction 


Both  polynomial  approximation  and  TPS  interpolation  are  implemented  to  correct  pincushion  distortion. 
The  difference  between  approximation  and  interpolation  lies  in  the  fact  that  approximation  methods  gener¬ 
ate  transformations  that  map  all  the  control  points  close  to  their  correspondence,  so  that  the  summation  of 
displacements  achieves  a  global  minimum;  whereas  interpolation  methods  produce  transformations  where 
all  the  control  points  can  be  mapped  to  their  correspondence  exactly  [Daview  and  Samuels  96]  (Appendix 
7.5). 


In  the  following  experiments,  control  points  are  chosen  to  be  equally  distributed  across  the  whole  image.  A 
grid  template  is  generated  as  a  gold  standard  for  comparison.  Fig.  8.  displays  the  corrected  results  by  both 
methods  with  different  numbers  of  control  points  chosen.  It  is  not  that  apparent  as  to  which  graph  shows  bet¬ 
ter  result.  Therefore,  we  need  to  use  some  quantitative  measurement  to  compare  the  results. 

Two  numbers  are  calculated  to  compare  the  quality  of  the  corrected  images:  the  error  rate  (e)  and  the  cross 
correlation  coefficient  (p). 

The  error  rate  e  is  defined  by  Eq.  (1).  Since  the  template  is  a  binary  image  with  the  brightness  of  grid  lines  as 
Os,  and  other  places  Is,  the  definition  of  the  error  rate  in  Eq.  (1)  is  the  same  as  the  mean  square  error  (MSE). 

e  -  sum  °f  misaligned  pixels  .. . 

total  number  of  pixels  '  ' 

The  cross  correlation  coefficient  p  is  defined  as  Eq.  (2). 
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=  cov(x,y)  =  g[(x-g[x])(y-£[y])] 
Jvar(x)var(y)  jE[x-E[xl\E[y-E[y]\ 


(2) 


Fig.  9  and  Fig.  10  are  two  plots  of  error  rate  comparison  and  cross  correlation  coefficients  comparison. 


Figure  8.  Correction  results  from  both  polynomial  approximation  (left  column)  and  TPS  interpolation  (right  column) 
with  respect  to  different  numbers  of  control  points.  From  top  to  bottom:  33  x  33, 17  x  17, 9  x  9,  and  5x5.  The 
artifacts  in  (b)  and  (d)  are  explained  on  page  9. 
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Figure  9.  Comparison  of  error  rates  of  polynomial  approximation  and  TPS  interpolation  as  a  function  of  the  number 
of  control  points. 


Figure  10.  Comparison  of  cross  correlation  coefficient  between  different  correction  results  and  the  template. 

From  the  above  comparisons,  we  draw  the  following  conclusions:  (1)  In  general,  TPS  interpolation  performs 
better  than  polynomial  approximation.  Except  that  with  5x5  control  points,  TPS  has  a  lower  error  rate  and 
a  higher  cross  correlation  coefficient  than  those  of  the  polynomial  approximation.  (2)  The  quality  of  the  cor¬ 
rected  image  relates  to  the  number  of  control  points  being  chosen.  Surprisingly,  it  is  not  true  that  the  more 
control  points,  the  less  the  error  rate.  The  corrected  image  has  the  lowest  error  rate  and  highest  correlation 
coefficient  with  9x9  control  points.  The  reason  is  stated  as  follows: 
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•  The  interpolation  approach  can  achieve  exact  mapping  based  on  the  assumption  that  the  control  point 
positions  are  known  exactly .  In  real  applications,  however,  the  positions  of  the  control  points  in  the  mea¬ 
sured  image  can  only  be  determined  approximately.  Therefore,  the  more  control  points  chosen,  the  more 
the  potential  error  might  accumulate,  resulting  in  the  severe  local  distortion  in  the  corrected  images  when 
choosing  33  x  33  control  points,  and  17  x  17  control  points,  as  can  be  observed  in  Fig.  8.  More  local  dis¬ 
tortions  exist  in  the  correct  image  with  33  x  33  control  points  than  that  with  17  x  17  control  points. 

•  For  approximation  approach,  the  more  control  points  chosen,  the  more  tendency  the  mapping  functions 
show  towards  discontinuity.  In  other  words,  the  mapping  function  might  fit  well  for  the  control  points, 
but  not  globally  for  all  the  points  in  the  image.  Therefore,  unless  we  choose  enough  control  points  (e.g. 
all  the  points  in  the  image),  or  we  will  get  more  errors  when  continue  increasing  the  number  of  control 
points. 

2.3  Missing  Data  Estimation 

Based  on  the  degradation  analysis  in  Appendix  7.6,  the  image  restoration  problem  can  be  formulated  as  Fig. 
11.  An  original  image /is  slightly  blurred  by  a  point  spread  function  (PSF)  h  from  the  X-ray  source.  The 
blurred  image  is  further  corrupted  by  fixed  noise  ( n )  and  defects  ( d)  from  the  CCD  detector,  g  is  the  so  called 
measured  image . 


original 

image 


blur 

® 


♦ 


h 


detector 


noise  induced 


n  d 


recovered 


missing  data  estimation  with 
deblurring  and  denoising 


measured 


image 


Figure  11.  System  model. 

The  image  formation  process  can  be  expressed  by  Eq.  (3),  where  0  stands  for  the  convolution  operator. 

g  =  f®h  +  n-d  (3) 

where  we  assume  that  the  PSF  is  already  known,  the  noise  distributes  as  independent  Gaussian,  and  the  posi¬ 
tions  of  missing  data  can  be  identified. 

We  propose  two  approaches  to  estimate  missing  data  along  with  deblurring  and  denoising:  the  consistency 
method  using  separable  deblurring,  and  a  maximum  a-posteriori  probability  (MAP)  method,  in  which  the 
optimization  is  solved  using  Mean  Field  Annealing  (MFA). 

The  basic  idea  behind  both  of  our  approaches  is  to  make  use  of  the  point  spread  function  (PSF):  before  a 
pixel  is  missed,  it  has  already  distributed  its  information  to  its  neighbors  through  the  effect  of  blur.  Fig.  12 
illustrates  this  process,  where  each  block  indicates  a  pixel,  and  the  blocks  marked  with  “0”s  are  the  missing 
pixels.  The  figure  shows  that  before  pixel  C  is  missed,  it  spreads  its  information  to  the  neighbors  gl,  g3,  g4, 
g6,  g7,  and  g9. 
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Figure  12.  Information  distribution  by  blur. 

If  we  can  obtain  a  measurement  of  the  PSF,  in  theory,  we  can  reconstruct  the  missing  data.  Regularization 
theory  (Appendix  7.7)  is  used  to  convert  an  ill-posed  problem  to  a  well-posed  one.  The  consistency  method 
using  separable  deblurring  (Appendix  7.9)  constructs  a  well-posed  problem  by  putting  several  restrictions  on 
the  image  data  and  the  blur  kernel.  The  MAP  method  (Appendix  7.10)  is  actually  a  Bayesian  interpretation 
of  regularization  problems,  which  is  solved  by  a  global  optimization  technique,  called  mean  field  annealing 
(MFA). 

Experiments  are  carried  on  both  synthetic  images  and  real  images.  System  distortions  and  degradations  are 
characterized  in  Appendix  7.8. 

The  synthetic  image  are  generated  by  simulating  the  image  formation  process:  (1)  blur  the  original  image  (/) 
by  convolving /with  a  finite  blur  kernel  h  calculated  in  Appendix  7.8;  (2)  insert  Gaussian  white  noise  with 
zero  mean  and  standard  deviation  equal  to  3.9,  which  is  obtained  in  Appendix  7.8;  and  (3)  remove  one  col¬ 
umn  of  data  to  generate  the  measured  image  g.  Five  synthetic  images  are  generated  showing  different 
degrees  of  smoothness  and  different  image  properties.  These  images  are  a  piecewise  constant  image,  a  piece- 
wise  linear  image,  a  piecewise  quadratic  image,  a  sinusoid  image,  and  a  mammogram  (Fig.  13). 


Figure  13.  Five  synthetic  images  with  different  properties.  From  left  to  right:  piecewise  constant,  piecewise  linear, 
piecewise  quadratic,  sinusoid,  and  mammogram. 

We  conduct  experiments  by  relaxing  the  following  two  assumptions  one  by  one  in  order  to  investigate  the 
effect  of  each  assumption  to  the  quality  of  the  restored  image: 

1)  the  blur  kernel  is  exactly  known;  and 

2)  noise  is  not  inserted. 

Experimental  results  from  both  the  consistency  method  (using  integer  criterion  and  neighbor  least  square 
error  criterion)  and  the  MFA  are  exhibited  and  compared  based  on  line  profiles  and  column  profiles. 

2.3.1  Results  Based  on  Complete  Set  of  Assumptions 

With  the  complete  set  of  assumptions  (the  blur  kernel  is  exactly  known,  the  blur  kernel  is  separable,  the  orig¬ 
inal  image  is  of  integer  type,  and  no  noise  is  inserted)  satisfied,  the  consistency  method  using  integer  crite¬ 
rion  can  recover  the  missing  column  exactly  from  the  measured  image.  The  consistency  method  using 
neighbor  least  square  error  (NLSE)  criterion  shows  better  results  than  the  MFA  method  in  deblurring,  as 
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Figure  15.  Profile  comparisons  of  restoration  results  from  the  five  synthetic  images.  In  each  plot,  from  top  to  bottom: 
profile  from  the  original  image,  restored  images  by  consistency  (int),  consistency  (NLSE),  and  MFA, 
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A  cross-correlation  comparison  between  the  original  image  and  the  restored  one  using  different  methods  is 
another  way  to  compare  the  algorithm  performance,  as  shown  in  Fig.  16. 


Figure  16.  Cross-correlation  comparison. 

The  consistency  method  with  integer  criterion  always  has  the  highest  cross-correlation  coefficient,  which  is 
1.  The  consistency  method  with  NLSE  criterion  shows  better  correlation  than  the  MFA  method  most  of  the 
time  except  for  sinusoid  images. 

2.3.2  Relaxing  the  Assumption  of  Perfect  Blur  Kernel  Estimation 

We  use  blur  kernel  h  to  blur  the  original  image,  but  use  h  +  8  to  restore  the  measured  image,  where  £  repre¬ 
sents  the  error  added  to  a  blur  kernel.  The  accuracy  of  an  estimated  blur  kernel  is  measured  by  correct  num¬ 
ber  of  digits.  For  example,  £  =  1%  means  the  2nd  digit  of  elements  in  the  blur  kernel  is  not  correct,  and  £  = 
0.1%  means  the  error  shows  up  in  the  3rd  digit.  Therefore,  if  h  =  [0.0625  0.1250  0.625  0.1250  0.0625],  a 
possible  estimation  when  £  =  1%  is  h  =  [0.07  0.13  0.63  0.13  0.07],  where  the  2nd  digit  is  not  accurate.  Fig. 
17  shows  restoration  results  of  consistency  method  using  NLSE  criterion  and  MFA  method  at  different  error 
rate  £  (0.1%  and  1%),  Fig.  18  is  a  column/row  profile  comparison,  and  Fig.  19  is  the  cross-correlation  com¬ 
parison. 
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Figure  17.  Restoration  results  with  an  inaccurate  blur  kernel 


Figure  18.  Profile  comparisons  of  restoration  results.  In  each  plot,  from  top  to  bottom:  profile  from  the  original 
image,  restored  images  by  consistency  (NLSE)  with  error  0.1%,  by  MFA  with  error  0.1%,  by  consistency  (NLSE) 
with  error  1%,  and  by  MFA  with  error  1%. 
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Figure  19.  Cross-correlation  comparison. 


From  Fig.  19,  we  can  see  that  the  MFA  method  is  not  sensitive  to  very  small  errors  in  the  blur  kernel  since 
the  correlation  coefficients  do  not  change  when  increasing  the  error  rate  from  0.1%  to  1%.  However,  the 
coefficients  of  the  consistency  method  with  NLSE  criterion  drops  when  increasing  the  error.  When  error  is 
on  the  2nd  digit  (e  =  0.1%),  the  consistency  method  provides  better  correlation  than  the  MFA  method  except 
for  the  sinusoid  image. 


To  summarize,  when  the  estimated  blur  kernel  is  in  high  accuracy  (two  or  more  than  two  digits  are  correct), 
the  consistency  method  performs  better  than  the  MFA  method.  However,  the  consistency  method  is  more 
sensitive  than  MFA  when  error  in  the  estimated  kernel  is  increased. 


2.3.3  Relaxing  the  Noise-Free  Assumption 

Gaussian  white  noises  nx  ~  N( 0, 0.012) ,  n2  ~  N( 0, 0.12)  and  n3  ~  N(0, 1)  are  inserted  to  the  blurred  image, 

and  one  column  of  data  are  deleted.  Restoration  results  on  the  piecewise  constant  image  are  shown  in  Fig. 
20.  Column/row  profiles  are  compared  in  Fig.  21  and  Fig.  22. 
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Figure  22.  Row  profiles  comparison. 

The  results  show  us  that  when  noise  is  very  small,  the  consistency  method  performs  better  than  the  MFA 
method,  it  achieves  a  higher  correlation  coefficient  (0.998269)  than  the  MFA  method  (0.998084).  However, 
when  increasing  the  standard  deviation  of  the  inserted  noise,  the  consistency  method  is  much  more  disturbed 
than  the  MFA  method. 
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2.3.4  Restoration  Results  of  Real  Image 


Fig.  23  exhibits  a  stripe  of  the  grid  template  image  (Fig.  44),  and  its  restoration  results  using  both  the  consis¬ 
tency  method  and  the  MFA  method. 


Figure  23.  Restoration  results  of  a  stripe  of  grid  template  image  with  flash  light  source.  From  left  to  right:  the 
measured  image,  restored  image  by  NLSE,  restored  image  by  MFA. 


The  restored  image  from  the  consistency  method  is  very  similar  to  that  from  the  MFA  method,  but  the  con¬ 
sistency  method  is  more  efficient  according  to  the  complexity  analysis  in  Appendix  7.1 1. 


To  summarize,  the  consistency  method  works  better  when  noise  is  very  small,  and  blur  kernel  is  estimated 
with  high  accuracy.  Performance  of  the  consistency  method  may  be  largely  affected  by  perturbations  in  the 
estimated  blur  matrix,  or  by  noise,  with  perturbations  in  noise  playing  a  more  important  role.  The  MFA 
method  behaves  more  stable  for  large  noise  or  inaccuracy  in  the  blur  kernel  estimation. 


2.4  Original  Contribution 

The  original  contribution  of  this  research  work  includes  the  use  of  conditioning  theory  in  the  analysis  of 
image  restoration  problem,  the  development  of  consistency  method  with  separable  deblurring,  and  the  use  of 
thin-plate  spline  interpolation  algorithm  in  geometric  correction. 

2.4.1  Conditioning  Analysis  of  Missing  Data  Estimation 

While  there  have  been  many  articles  in  the  literature  discussing  various  image  restoration  algorithms,  few 
[Forbes  and  Anh  94][Milinazzo  et  al.  87]  pay  attention  to  quantify  the  ill-conditioning  of  the  blur  kernel. 
The  most  important  contribution  of  this  research  work  is  that  we  use  conditioning  theory  to  analyze  image 


July  15,  1999 


19 


restoration  problems,  which  has  helped  us  obtain  an  in-depth  view  of  how  algorithms  might  behave  theoret¬ 
ically  and  how  perturbation  of  each  parameter  might  affect  the  solution.  We  also  use  the  conditioning  analy¬ 
sis  to  guide  our  algorithm  design  and  parameter  selection.  The  conclusions  we  draw  from  the  conditioning 
analysis  are  well  demonstrated  by  the  experimental  results. 

2.4.2  Consistency  Method  with  Separable  Deblurring 

Neither  separable  convolution  nor  the  idea  of  converting  a  convolution  operation  into  a  matrix  multiplication 
is  new.  These  techniques  have  been  addressed  very  early,  back  to  70s  [Andrews  and  Hunt  77].  Our  contribu¬ 
tion  is  that  we  first  take  advantage  of  the  separability  property  of  some  blur  kernel  (like  Gaussian),  and  con¬ 
vert  the  separable  convolution  operation  into  a  separable  matrix  multiplication.  By  inserting  several 
restrictions  on  the  data,  we  can  construct  the  separated  circulant  matrix  and  be  able  to  demonstrate  that  the 
condition  number  of  such  matrices  is  very  small,  so  that  the  problem  is  well-conditioned.  In  this  way,  we 
convert  the  ill-posed  image  restoration  problem  to  a  well-posed  one. 

By  using  a  backward-stable  algorithm,  we  are  able  to  restore  the  missing  data  with  100%  accuracy  assuming 
the  defects  only  come  from  blur  and  the  original  image  is  of  integer  type. 

The  beauty  of  this  algorithm  lies  in  its  simplicity  and  computational  efficiency. 

2.4.3  TPS  in  Geometric  Correction 

TPS  interpolation  is  a  conventional  method  for  interpolation.  We  first  apply  it  to  correct  the  geometric  dis¬ 
tortion.  No  other  published  literature  has  used  this  method  under  this  scenario.  By  using  TPS  interpolation, 
we  gain  better  performance  in  the  correction  results  than  the  traditional  polynomial  method. 
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7.  Appendixes 

7.1  Image  Acquisition 

The  image  acquisition  module  consists  of  three  parts:  (1 )  the  CCD  head  module ,  that  converts  optical  signals 
to  electrical  signals;  (2)  the  CCD  driver  assembly ,  that  provides  the  system  control  signals,  the  sampling 
mechanisms,  and  the  signal  amplification  circuits;  and  (3)  the  analog-to-digital  converter  (ADC),  that  quan¬ 
tizes  and  resamples  (optional)  analog  signals  and  converts  them  to  digital  ones. 


Figure  24.  Front  view  of  the  image  acquisition  system  assembly. 
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Figure  25.  System  connection  diagram. 

Fig.  24  shows  a  front  view  of  the  image  acquisition  system  assembled  in  a  Bud  box.  Rash  is  used  as  the  light 
source  instead  of  X-ray  for  experimental  purpose.  Fig.  25  is  a  simplified  system  hookup  diagram. 

7.1.1  CCD  Head  Module 

The  head  module  (from  EEV,  Inc.)  performs  like  a  detector  and  a  transducer.  It  contains  a  scintillator  screen, 
a  1  x  2  CCD  sensor  arrays,  each  of  which  is  hard  coupled  to  a  fiber  optic  taper.  The  two  tapers  are  dry  cou- 
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pled  together  in  a  mechanical  frame,  so  that  the  dead  space  between  adjacent  tapers  is  always  less  than 
200|im,  with  50pm  the  usual  space. 

missing  information  about  50  micron  width 


Figure  26.  The  CCD  head  module. 

The  scintillator  uses  a  phosphor  screen  to  convert  an  X-ray  beam  into  myriad  points  of  light.  The  phosphor 
screen  is  Eu  doped  Gd202S,  100pm  thick. 

Each  fiber  optic  taper  has  a  50  x  50  mm2  outer  surface  and  a  25  x  25  mm2  inner  surface,  which  is  usually  less 
than  the  CCD  active  area.  The  CCD  active  area  is  25.9  x  27.9  mm2,  consists  of  1152  rows  and  1242  col¬ 
umns,  with  each  pixel  size  22.5  x  22.5  pm2. 

Light  is  guided  from  the  bottom  of  phosphor  screen  down  to  each  CCD  cell  by  a  bundle  of  fibers.  The  CCD 
cell  can  then  convert  the  light  signal  to  an  electrical  signal  by  photoelectric  process,  and  the  generated 
charges  are  accumulated  at  each  CCD  cell. 

Care  must  be  taken  not  to  allow  too  much  light  onto  the  CCD  since  a  gross  overload  can  cause  the  resulting 
image  to  turn  black.  Basically,  too  much  light  entering  the  CCD  causes  the  video  signal  out  of  the  device  to 
look  like  there  is  no  light  entering  the  CCD  (i.e.  black).  However,  with  no  light  entering  the  front  of  the  optic 
taper  but  with  light  still  able  to  enter  through  the  sides,  there  is  still  enough  light  entering  the  CCD  to  cause 
the  video  signal  out  of  the  CCD  to  become  saturated  (i.e.  white).  Thus  careful  control  of  the  lighting  condi¬ 
tions  is  necessary  to  achieve  correct  results.  In  our  experiments,  the  Bud  box  is  covered  by  opaque  black 
paper  and  the  edges  of  the  fiber  optic  tapers  are  covered  by  black  insulating  tape  to  prevent  light  leakage  into 
the  CCD.  The  flash  light  source  is  covered  by  two  layers  of  black  insulating  tapes,  with  a  0.5  x  0.5  cm2 
square  hole  in  the  outer  layer  tape  as  our  point  source  to  reduce  the  amount  of  light  getting  to  the  CCD. 

7.1.2  CCD  Driver  Assembly 

CCD  driver  assembly  (from  EEV,  Inc.)  drives  the  CCD  readout.  Three  boards  are  involved  in  this  process: 
digital  board,  analog  board,  and  buffer  board  (Fig.  25). 

Digital  Board.  The  digital  board  houses  the  control  logic  for  the  CCD  sensor  and  the  video  processing.  It  is 
also  responsible  for  interfacing  with  the  outside  world.  The  digital  board  sends  control  pulses  to  the  analog 
board  via  a  20-way  piggyback  connector,  and  to  the  buffer  board  via  a  two-meter  37-way  ribbon  cable.  The 
64- way  (a  and  c)  DIN  41612  connector  is  where  power  is  fed  into  the  CCD  Driver  Assembly  and  where  the 
control  and  synchronization  signals  enter  and  leave.  The  DIN  connector  we  use  is  from  Harting  (0903  164 
6821).  It  is  a  female  connector  with  13mm  wrap  posts,  64  contacts  mounted  on  columns  a  and  c. 

Four  of  the  most  important  control  signals  are  pixel  clock ,  line  blank ,  grab ,  and  request  integration  (RI), 
where  the  line  blank  serves  as  the  horizontal  synchronization  signal,  and  the  grab  as  the  vertical  synchroni¬ 
zation  signal.  The  timing  diagram  is  shown  in  Fig.  27. 
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•  Line  blank  is  a  pulse  that  occurs  once  every  line  read  out  from  the  CCD.  It  is  LOW  during  the  output  of 
active  video,  HIGH  during  line  transfer,  readout  of  blank  elements,  dark  reference  elements  and  transi¬ 
tion  elements  (Fig.  27  (a)). 

•  Grab  is  used  to  identify  active  video  from  the  CCD.  Grab  is  LOW  during  the  readout  of  active  lines, 
HIGH  at  all  other  times  (Fig.  27  (b)). 

•  Request  for  integration  (RI)  determines  the  integration  interval.  When  RI  goes  HIGH,  clocking  stops, 
integration  starts;  when  it  goes  LOW,  readout  starts  (Fig.  27  (b)).  RI  triggers  the  Grab  signal.  Fig.  28  is 
the  circuit  design  for  generating  RI  signal.  The  integration  interval  is  about  100  ms. 

Buffer  Board  and  Analog  Board.  The  buffer  board  terminates  the  control  signals  from  the  digital  board  to 
the  sensor  and  amplifies  the  output  signal  from  the  sensor.  The  sensor  signal  is  sent  via  a  mini-COAX  cable 
(SMC  port)  to  the  analog  board  (Fig.  24,  Fig.  25). 

The  analog  board  processes  the  sensor  signal  from  the  buffer  board  using  correlated  double  sampling  (CDS) 
to  isolate  the  signal  level  difference,  and  produce  analog  video,  with  adjustable  gain.  The  video  signal  is 
linked  to  a  BNC  connector  on  the  front  panel  which  will  go  into  the  input  channel  of  the  A/D  board  from 
GaGe. 

Amplification  and  Readout  Settings.  The  signal  read  out  from  the  CCD  head  module  is  very  weak,  and 
amplification  must  be  conducted  before  digitization.  In  this  driver  assembly,  gains  are  set  at  three  stages  of 
the  video  path.  The  first  gain  occurs  on  the  buffer  board,  where  an  amplifier  functions  with  adjustable  gains 
at  4,  8,  or  16.  The  second  is  on  the  analog  board  with  adjustable  gains  at  1  or  4  that  adjusts  the  input  to  the 
CDS  circuit.  The  last  gain  stage  is  under  control  of  a  parameter  switch  on  the  digital  board,  that  adjusts  the 
output  from  the  CDS  circuit  at  xl,  x2,  or  x4. 

There  are  6  switches  on  the  digital  board  (Fig.  29),  three  of  which  are  used  to  set  up  important  parameters 
that  control  the  CCD  readout.  Following  is  an  overview  of  these  three  switches,  and  the  appropriate  settings 
for  our  system. 
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(on  front  panel,  reset  button) 
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Switch  1 
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Figure  29.  Switches  on  the  digital  board. 


Switch  1  controls  the  readout  mode,  and  the  readout  rate.  Switch  2  sets  up  the  gains  for  final  video  readout. 
Switch  5  is  a  front  panel,  reset  button,  which  resets  the  power  supply  of  the  driver  board.  Switch  6  indicates 
the  model  of  the  CCD  assembly.  Table  2  is  the  settings  for  each  switch,  and  Table  3  is  a  brief  explanation  of 
the  parameters  that  the  switch  controls.  More  details  are  addressed  in  [EEV  94]. 
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TABLE  2.  Switch  settings. 


TABLE  3.  Explanation  of  switch  parameters. 


■ 

Switch  1 

Switch  2 

on 

off 

on  off 

a 

Frame  transfer 

Full  frame 

Direction  of  vertical  readout  charge  (0: 
down;  1:  up) 

n 

Standard  mode 

Inverted  mode 

c 

Readout  rate 

The  setting  for  1MHz  is  010 

Direction  of  horizontal  readout  charge  (01: 
right;  10:  left) 

d 

e 

Gain  setup  for  output  of  CDS  circuit  (00: 
x8,  10:  x4,  01:  x2,  00:  xl) 

f 

Normal 

TDI  mode 

g 

Remote/Computer 

Local  Switch 

N/A 

mm 

Visible  light 

X-ray  source 

N/A 

7.1.3  Other  Supply  Accessories 

Besides  the  CCD  driver  assembly,  other  supply  accessories  are  needed  for  successful  CCD  readout.  These 
accessories  include  RI  circuit  board,  power  supply,  test  pod,  etc. 

Request  integration  (RI)  is  a  key  control  signal  in  the  driver  assembly.  When  RI  is  high,  the  CCD  sensor 
stops  clocking,  allowing  light  accumulation,  which  is  the  integration  process.  After  a  100ms  integration 
interval,  RI  goes  LOW,  and  drives  GRAB  to  go  HIGH  to  start  the  readout  process.  RI  circuit  board  (Fig.  28) 
is  designed  to  generate  a  100ms  pulse,  controlled  by  an  SPDT  momentary  switch. 

To  synchronize  the  start  of  integration  and  the  flash,  the  RI  circuit  board  adopts  some  extra  design  ideas  to 
achieve  high-speed,  noise-free  logic  [Snyder  93].  A  large  tantalum  capacitor  (10  to  20  pF)  bypasses  the 
board  to  suppress  the  low-frequency  external  noise.  Each  chip  on  board  is  bypassed  with  0.001  to  0.01  pF 
ceramic  or  monolithic  capacitors  between  Vcc  and  ground  to  avoid  sudden  voltage  drop  when  chips  demand 
sudden  bursts  of  current. 

To  support  the  CCD  head  module  and  driver  assembly,  a  power  supply  is  required  to  provide  +5  volts  at 
400mA,  +15  volts  at  500mA,  and  -15  volts  at  300mA.  For  a  1  x  2  CCD  sensor  arrays,  in  order  to  drive  both 
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of  the  CCD  sensors  at  the  same  time,  the  current  requirements  for  each  voltage  level  are  doubled.  That  is,  +5 
volts  at  800mA,  +15  volts  at  1  A,  -15  volts  at  600mA. 


Figure  30.  Power  supply  (International  PowerPC  Power  Supplies  -  IHBCC-512)  connection  diagram. 

The  power  supply  we  used  is  from  International  Power  DC  Power  Supplies.  The  IHBCC-512  model  pro¬ 
vides  three  outputs:  5V  @  3A,  15V  @  3 A,  and  -15  V  @  3 A,  which  is  enough  for  our  application.  Fig.  30 
shows  the  connection  diagram  inside  the  power  supply. 


July  15, 1999 


35 


Test  pod  (Fig.  31)  is  used  to  test  the  signals  from  the  DIN  connector  that  are  generated  by  the  digital  assem¬ 
bly.  It  is  also  used  to  switch  between  channel  1  control  signals  and  channel  2  control  signals  to  correctly  trig¬ 
ger  the  ADC. 


from  DIN  connector 


Figure  31.  Test  pod  signals. 

7.1.4  Analog-to-Digital  Converter 

The  GaGe  data  acquisition  board  functions  as  a  12  bit  A/D  converter  which  can  quantize  and  resample  the 
analog  signal  into  a  12-bit  digital  signal.  The  maximum  data  transfer  rate  is  approximately  1.5  MBytes/sec. 


Trigger  Address 


Figure  32.  The  circular  on-board  memory  of  the  Gage  data  acquisition  board. 
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The  data  coining  out  of  the  ADC  is  stored  in  the  static  memory  (SRAM)  of  the  GaGe  board  which  is  config¬ 
ured  as  a  circular  buffer  (Fig.  32  (a)).  Circular  buffer  is  used  to  guarantee  that  the  system  will  keep  on  cap¬ 
turing  data  indefinitely  until  a  trigger  event  is  detected. 

Once  initialized,  the  board  starts  to  read  in  data,  and  at  the  same  time  waiting  for  trigger  event  to  happen. 
The  data  that  read  in  during  this  period  is  called  pre-trigger  data.  When  trigger  event  is  received,  a  specified 
number  of  pixels  are  captured.  These  data  are  called  post-trigger  data.  After  storing  the  specified  number  of 
post-trigger  data,  the  acquisition  is  stopped  (Fig.  32  (b)). 

The  board  we  used  is  model  CSS  12,  which  must  be  working  under  dual  channel  mode.  That  is,  the  two  chan¬ 
nels  of  data  can  be  input  at  the  same  time.  The  array  elements  location[ 0]  and  location[  1]  hold  the  first  pixel 
value  from  channel  1,  with  location[  1]  storing  the  least  significant  8  bits,  and  location[0 ]  storing  the  most 
significant  8  bits.  Similarly,  location[ 2]  and  location [3]  hold  the  first  pixel  value  from  channel  2. 

The  sample  rate  of  the  ADC  is  set  at  2MHz,  with  each  channel  at  1MHz,  thus  matching  the  CCD  readout 
rate.  Since  the  on-board  memory  of  this  model  is  1M,  with  512k  for  each  channel,  three  integrations  and 
readouts  are  required  to  capture  images  from  one  CCD  sensor  array,  and  six  from  both  arrays. 


With  a  1MHz  readout  rate,  1M  sample  rate,  1152  vertical  pixels  and  1242  horizontal  pixels  (equating  to 
1,430,784  total  pixels  in  an  array),  the  performance  of  our  data  acquisition  module  can  be  summarized  as 
follows:  1.43s  minimum  is  required  to  readout  a  complete  frame;  1.24ms  minimum  is  required  to  readout  a 
complete  line;  and  l|is  is  required  to  readout  1  pixel  of  information. 

7.2  Image  Display 

Image  display  module  may  be  needed  anywhere,  from  measured  image  display,  corrected  image  display,  to 
restored  image  display. 

In  our  system,  the  spatial  resolution  of  a  complete  image  is  roughly  1152  x  2484,  involving  a  1  x  2  1152  x 
1242  arrays.  Since  most  of  the  current  commercial  monitors  can  only  display  image  with  dimensions  1024  x 
1024  or  smaller,  how  to  efficiently  and  conveniently  display  larger  images  on  smaller  display  planes  is 
another  problem  which  needs  to  be  solved. 

Here,  we  first  discuss  two  possible  approaches  for  image  display:  (1)  use  the  on-board  video  card ;  and  (2) 
use  a  secondary  display  board.  We  then  propose  to  combine  these  two  approaches  and  achieve  dual-screen 
mode  loupe  display  in  our  system. 

7.2.1  Image  Display  through  On-Board  Video  Card 

This  is  the  most  common  approach  for  image  display.  Since  our  A/D  board  works  only  under  Windows 
operating  system,  we  choose  Visual  C++  as  the  programming  language,  and  use  VisSDK  (Vision  Software 
Developers  Kit)  as  the  display  library.  VisSDK  is  a  low-level  C++  library,  developed  by  Microsoft  Vision 
Group.  It  aims  to  provide  a  strong  programming  foundation  for  research  and  application  development,  and 
includes  possible  real-time  imaging.  VisSDK  is  organized  into  different  projects  and  VisDisplay  is  the  one 
covers  display  related  functions. 

There  are  two  ways  to  display  a  single  image  in  VisDisplay.  We  can  either  get  access  to  the  raw  pixel  data 
and  build  a  Windows  bitmap  from  it,  or  we  can  use  a  Windows  HDC  (handle  of  a  device  context).  A  device 
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context  (DC)  is  essentially  a  unified  interface  to  a  specific  graphical  device.  It  can  be  used  to  obtain  informa¬ 
tion  about  the  device  and  also  perform  graphical  output  to  the  device. 


Figure  33.  Software  layers  in  image  display. 

Fig.  33  is  a  software  layer  diagram.  At  the  lower  level,  device-dependent  drivers  are  used  to  operate  on  the 
corresponding  hardware  directly.  At  a  higher  level,  the  graphics  device  interface  (GDI)  performs  primitive 
device-independent  graphical  functions  through  DC.  MFC  (Microsoft  Fundamental  Class)  is  another  layer 
that  is  built  upon  GDI.  It  consists  of  a  series  of  predefined  classes  which  save  the  users  from  tedious  primi¬ 
tive  GDI  function  calls,  and  provide  a  more  convenient  interface. 

Once  the  image  is  read  into  the  host  memory,  it  can  be  displayed  in  a  window’s  client  area.  Scrolling  in  both 
horizontal  and  vertical  directions  is  needed  for  large  image  display.  A  window’s  client  area  is  a  versatile  sur¬ 
face  that  can  display  anything  a  Windows  program  can  draw. 

7.2.2  Image  Display  by  Secondary  Display  Board 

The  secondary  display  board  we  used  is  the  Matrox  Pulsar  from  Matrox  Electronic  Systems  Ltd.  This  board 
features  both  on-board  grab  and  display  capabilities.  The  reason  we  did  not  use  its  grab  functionality  is 
because  the  minimum  pixel  clock  it  can  accept  is  2MHz,  while  the  pixel  clock  from  our  image  acquisition 
system  is  1MHz.  However,  the  Matrox  board  provides  a  powerful  display  section,  featuring  a  2M  image 
frame  buffer,  a  2M  graphics  overlay  (VGA)  frame  buffer  for  non-destructive  overlay  capabilities,  and  dual¬ 
screen  mode  display. 
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To 

Monitor 


Secondary  32-bit  PCI  local  BUS 
Figure  34.  Display  section  of  Matrox  Pulsar. 

The  display  section  (Fig.  34)  is  also  powered  by  the  Matrox  MGA-2064W  graphics  accelerator  that  can  dis¬ 
play  at  resolutions  up  to  1600x1200x256  colors.  This  graphics  accelerator  can  be  used  either  as  the  main  dis¬ 
play  controller  of  the  host  system;  or  as  a  separate  graphics  board  so  that  users  can  work  in  dual-screen 
mode.  The  overlay  (VGA)  and  main  image  frame  buffers  support  common  zooming  (by  factors  of  2  or  4), 
panning,  and  scrolling. 

A  software  package  called  MIL  (Matrox  Imaging  Library),  developed  by  Matrox,  provides  interfacing  func¬ 
tion  calls  to  the  display  board. 

The  big  advantage  of  using  a  second  display  board  is  the  dual-screen  display  mode.  The  users  can  display 
images  on  one  monitor,  while  doing  controls  on  another. 


7.2.3  Dual-Screen  Mode  Loupe  Display 

In  medical  imaging,  very  often  some  detail  analysis  is  needed  on  a  few  interesting  sections  of  the  image, 
such  as  a  suspicious  lesion  area.  For  the  purpose  of  conveniently  processing  sections  of  images,  while  at  the 
same  time  keeping  a  complete  viewgraph  of  the  original  one,  we  propose  to  combine  the  use  of  both 
approaches  we  discussed  above,  and  implement  the  so-called  dual-screen  model  loupe  display. 

We  use  the  word  loupe  because  it  is  analogous  to  the  operation  of  zooming  a  patch  of  the  original  image.  To 
keep  viewgraphs  of  images  at  different  scales,  dual-screen  is  a  natural  choice.  The  Matrox  Pulsar  board  also 
provides  built-in  functions  for  zooming  and  color  mapping,  that  can  be  easily  used  to  operate  on  the  zoomed 
patches  of  the  original  image. 

Object-oriented  programming  language  is  used  in  developing  this  technique.  Three  objects  are  extracted:  the 
Gage  ADC,  the  Pulsar  display  board,  and  the  image .  Gage  object  consists  of  methods  for  acquiring  data 
from  the  CCD  sensor  arrays.  Image  object  involves  methods  for  image  display  on  the  host  monitor,  and  for 
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image  segmentation.  Pulsar  object  is  responsible  for  displaying  segments  of  images  onto  the  secondary 
monitor,  along  with  methods  for  zooming,  color  mapping,  and  contrast  tuning,  etc.  (Fig.  35). 


Figure  35.  Object-oriented  display  design. 


We  developed  our  display  system  on  Windows  NT  platform.  Images  are  captured  using  the  interfacing  func¬ 
tions  provided  by  GaGe  ADC.  The  complete  image  is  displayed  on  the  main  monitor  (host  display).  The  sec¬ 
ondary  monitor  is  used  to  display  regions  of  interest  based  upon  user’s  selection.  Users  can  either  mark  the 
region  of  interest  by  a  rectangular,  or  they  can  click  on  any  interesting  pixel  with  a  surrounding  region  of  64 
x  64  or  128  x  128  automatically  extracted  from  the  original  image. 

The  region  of  interest  is  displayed  on  the  secondary  monitor,  where  several  image  enhancement  techniques 
can  be  operated.  These  techniques  include  zooming  at  any  integer  scale,  autoscaling,  contrast  stretching, 
pseudo-color  enhancement  (use  color  to  help  visual  discrimination  of  small  grey  level  differences),  edge 
enhancement  (convolve  image  with  high-pass  filter  kernel  such  as  Laplacian  type  kernels),  unsharp  masking 
(equivalent  to  high-pass  spatial  filtering),  etc. 

7.3  Sources  of  Radiometric  and  Geometric  Distortions 

In  our  system,  the  functions  of  the  fiber-optic  tapers  are  three-fold:  (1)  used  as  an  aligner  to  make  adjacent 
CCD  detectors  aligned  very  close  to  each  other  such  that  only  1~2  columns  of  information  is  missed  along 
the  butting  edge;  (2)  used  as  the  light-guide  system,  where  a  bundle  of  fibers  connect  each  point  on  the  scin¬ 
tillator  screen  to  a  corresponding  cell  on  the  CCD  detector;  and  (3)  used  as  a  demagnifier  to  increase  the 
field  of  view  of  the  CCDs,  which  is  usually  in  the  ratio  of  2:1. 

However,  the  use  of  the  fiber-optic  tapers  also  causes  radiometric  and  geometric  distortions,  which  origi¬ 
nates  during  manufacture.  Fig.  36  illustrates  how  a  fiber-optic  taper  is  constructed.  A  fiber-optic  tube  is  first 
heated  from  the  middle;  then  compressed  to  the  middle,  such  that  the  area  of  the  cutting  surface  matches  the 
size  of  a  CCD  array  to  which  it  will  attach;  finally  a  cut  from  the  middle  makes  two  fiber-optic  tapers  with 
the  smaller  end  attaches  to  the  CCD  array  and  the  other  end  the  outer-surface,  facing  the  scintillator. 


Heat  from  the  middle 


Q^=) 

ozrO 

Cut  in  the  middle 


Figure  36.  The  manufacturing  of  the  fiber-optic  tapers. 

In  this  process,  we  can  see  that  no  matter  how  precise  the  compression  process,  the  cutting  surface  of  the 
taper  can  not  perfectly  match  the  size  of  the  CCD  array.  Usually,  the  taper  is  more  compressed  than  it  should 
be  in  order  to  avoid  missing  data  sensed  at  the  boundary  of  the  outer-surface  of  the  taper.  That  is  why  the 
captured  image  is  always  less  than  the  CCD  active  area.  The  compression  thus  causes  the  “pincushion  dis¬ 
tortion”,  a  well-known  type  of  geometric  distortion,  where  both  the  horizontal  and  vertical  lines  bend 
inwards  toward  the  center  of  the  display  (Fig.  38).  Instead  of  each  bundle  of  fiber  being  guided  to  a  corre- 
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sponding  cell  of  the  CCD  detector,  the  bundles  at  the  boundaries  of  the  taper  is  deformed  like  Fig.  37.  The 
distortion  is  not  symmetric,  exhibiting  some  degree  of  shear  effect. 


Figure  37.  The  correct  alignment  (left),  the  misalignment  of  fiber  optic  bundles  and  CCD  cells  along  the  edge  of  the 
detector  which  causes  the  pincushion  distortion  (right). 


Figure  38.  Vignetting-type  radiometric  distortion  and  geometric  distortion. 


Furthermore,  because  of  the  shape  of  the  taper,  the  fibers  at  the  boundary  need  to  travel  a  longer  distance 
than  those  at  the  center,  and  the  light  transfer  is  less  efficient  when  the  fiber  is  not  normal  to  the  CCD  sur¬ 
face.  Thus,  there  is  an  intensity  gradation  from  the  center  to  margin  like  a  cross  section  of  a  tree,  or  a  set  of 
concentric  ellipses  (Fig.  38).  This  kind  of  distortion  is  called  the  vignetting-type  radiometric  distortion,  as  it 
resembles  the  vignetting  which  results  from  an  imperfect  lens. 


7.4  Radiometric  Correction 

Radiometry  is  the  study  of  image  formation.  It  concerns  two  sources  of  light  energy  [Trucco  and  Verri  98]: 
the  amount  of  light  reflected  from  the  object,  and  the  amount  of  reflected  light  that  actually  reaches  the  sen¬ 
sor.  Assume  the  amount  of  light  reflected  from  the  object  is  homogeneous;  if  there  is  no  radiometric  distor¬ 
tion,  the  captured  image  should  have  homogeneous  illumination  over  the  whole  image.  However,  as 
mentioned  in  Appendix  7.3,  the  vignetting-type  radiometric  distortion  corrupts  the  homogeneity. 

The  fundamental  equation  of  radiometric  distortion  is  Eq.  (4)  [Trucco  and  Verri  98],  where  g(;c,  y)  and  r|) 
are  the  image  irradiance  (the  power  of  the  light  at  each  point  of  the  image  plane  per  unit  area)  and  the  image 
radiance  (the  power  of  the  light  ideally  emitted  by  each  point  of  a  surface  in  3-D  space  in  a  given  direction) 
respectively.  Eq.  (4)  says  that  the  image  irradiance,  g(x ,  y),  decreases  as  the  fourth  power  of  the  cosine  of  the 
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angle  (a)  formed  by  the  principal  axis  with  the  optical  image  irradiance  (Fig.  39).  g(jt,  y)  is  also  regarded  as 
uniformly  proportional  to  the  scene  radiance,  ffe,  rj),  over  the  whole  image  plane. 


Figure  39.  Radiometry  in  image  formation. 


4  Kf. 


cos4a~  /(4,Ti)cos4a 


(4) 


Conventional  method  [Frei  83]  simplifies  the  radiometric  model  of  Eq.  (4)  by  a  polynomial  approximation  as 
Eq.  (5), 

*(*.?)« — (5) 

a0  +  a j(x  +  y  ) 


We  observe  from  Fig.  38  that  the  illumination  distribution  of  both  the  left  and  right  part  of  the  image  shaped 
like  ellipses  with  different  length  of  axes.  Therefore,  we  modify  Eq.  (5)  into  Eq.  (6),  where  the  equation  for 
an  ellipse  is  used  to  model  the  distortion,  (jc,-,  yj)  are  the  coordinates  of  pixel  i,  and  ( x{)>  ya)  denotes  the  origin 
of  the  concentric  ellipses. 


g(*i>  J/) 


2  2 
a0  +  ax(Xi-xo)  +ay(.yi~y0) 


(6) 


The  coefficients  a  =  (a0>  ax-  ay)  niust  be  determined  in  order  to  be  used  to  correct  the  distorted  image  g  back 
into  /  (Eq.  (7)),  such  that  Eq.  (8)  can  reach  its  minimum.  M  is  the  total  number  of  pixels  in  the  image. 


h  =  g(*r  y,)[«o + ax(xi  -  x,f  +  ay(yi  -  y,/} 


M  - 


e  =  mina  S 
1  =  0 

To  make  the  computation  easier,  we  rewrite  Eq.  (7)  and  Eq.  (8)  using  matrix  notation.  Let 
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then  Eq.  (7)  can  be  written  as 

F  =  GWA 


and  Eq.  (8)  as  Eq.  (9), 


e  =  minA{F-GWA)T(F-GWA) 


The  solution  to  Eq.  (9)  is  Eq.  (10), 

A  =  (GW0~‘f 

where  (GW)-1  is  the  pseudo-inverse,  which  can  be  computed  by 

(GWO-1  =  [{GW)T(GW)]-'(GW)T 


(9) 

(10) 


In  real  world  applications,  we  do  not  solve  Eq.  (9)  based  on  the  whole  set  of  image  pixels,  which  would  be 
very  time  consuming  and  unnecessary.  Instead,  we  choose  a  set  of  m  pixels  as  the  control  points.  The  control 
points  are  the  ones  that  we  know  what  their  original  intensity  values  should  be. 

We  choose  a  rectangular  grid  as  the  testing  image  (Fig.  38).  If  not  for  the  radiometric  distortion,  blur  or 
noise,  the  measured  image  would  be  binary,  with  only  two  intensity  levels:  255  and  0.  We  neglect  noise  for 
the  time  being,  and  choose  the  center  point  of  each  square  as  the  control  points  to  avoid  blur  effects.  The 
original  intensity  value  of  these  control  points  should  be  equal  to  the  intensity  value  of  the  origin  of  the 
ellipse. 

With  the  control  points  chosen,  we  can  easily  construct  matrices  W>  G,  and  F,  and  apply  them  to  Eq.  (10)  to 
solve  the  coefficient  matrix  A.  With  the  coefficients  solved,  a  look-up  table  can  be  generated  with  each  entry 
indicating  the  scale  that  should  be  imposed  onto  each  pixel  of  the  measured  image  g(x,  y)  to  correct  the  radi¬ 
ometric  distortion. 

7.5  Geometric  Correction 

In  geometric  correction,  a  transformation  function  is  designed  to  map  the  control  points  from  a  known  pat¬ 
tern  (such  as  a  regular  grid)  to  their  measured  positions.  The  geometric  distortion  in  any  captured  images  can 
then  be  corrected  by  back-projection.  Control  points  are  points  whose  positions  in  the  original  scene  are 
known  a-priori.  Polynomial  approximation,  usually  in  the  2nd  degree,  is  the  conventional  transformation 
function  most  often  used  [Pratt  78][Rosenfeld  and  Kak  82]. 

In  our  imaging  system,  we  implemented  the  polynomial  method  in  a  higher  degree.  We  also  adopted  the 
thin-plate  spline  (TPS)  as  another  type  of  transformation  function.  TPS  is  a  classical  method  for  interpola- 
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tion.  Bookstein  [Bookstein  89]  first  used  TPS  to  model  the  deformations  between  two  sets  of  landmark 
points.  As  far  as  we  know,  this  thesis  is  the  first  application  of  the  method  to  geometric  correction. 


Figure  40.  The  transformation  system  between  the  measured  image  and  the  corrected  image. 

The  difference  between  approximation  and  interpolation  lies  in  the  fact  that  approximation  methods  gener¬ 
ate  transformations  that  map  all  the  control  points  close  to  their  correspondence,  so  that  the  summation  of 
displacements  achieves  a  global  minimum;  whereas  interpolation  methods  produce  transformations  where 
all  the  control  points  can  be  mapped  to  their  correspondence  exactly  [Daview  and  Samuels  96]. 


7.5.1  Polynomial  Approximation 

Nearly  all  the  methods  used  to  correct  pincushion  distortion  so  far  are  variants  of  polynomial  approximation 
[Beier  et  al.  92][Butler  and  Pierson  91].  Pratt  addressed  the  2nd  degree  polynomial  approximation  in  1978 
[Pratt  78].  Here,  we  generalize  this  method  to  an  nth  degree  polynomial. 


Assume  yj| T  =  [Px(uit  vf)  Py(uit  v()]  ? ,  where  (uh  v,)  is  a  control  point  from  the  corrected  image,  (*f ,  y{ ) 

is  the  corresponding  point  in  the  measured  image,  and  ( Px ,  Py)  are  the  transformation  functions  for  the  x  and 
y  coordinates  respectively,  both  of  which  are  nth  degree  polynomials  with  respect  to  u  and  v.  Px  and  Py  are 
expressed  in  Eq.  (11)  and  Eq.  (12), 

n-  1 

X  =  Px(u,v)=  Yj  X  air/vS  (U) 

/  s  0r  +  A’  =  * 

y  =  Py(u,v)  =  £  X  hirs**  (12) 

i=0r+s=i 

where  airs  and  birs  are  the  coefficients.  The  transformation  functions  should  map  all  control  points  (n,  v)  to 
(x,  y)  as  closely  as  possible,  that  is,  (*,  y)  as  close  to  (;t,  y)  as  possible.  Minimum  mean  square  error 
(MMSE)  is  used  to  measure  the  error  as  Eq.  (13), 

m  -  1 

E  =  mina  h  £  [(xrxi)2  +  {yryi)2]  (13) 

i  =  0 

where  m  is  the  number  of  control  points.  Coefficients  a  and  b  are  determined  so  that  £  can  reach  its  mini¬ 
mum. 


Since  the  pincushion  distortion  that  our  imaging  system  suffers  is  not  symmetric,  we  use  a  3rd  degree  poly¬ 
nomial  to  model  the  transformation.  Based  on  Eq.  (11)  and  Eq.  (12),  we  get 


2  2  3  2  2  3 

X  =  a0oo  +  amu  +  amv  +  a220u  +  a2nuv  +  a202v  +  amu  +a32lu  v  +  a3l2uv  +  amv 
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and 


y  =  bOOO  +  bmu  +  bmv  +  b220“  +  b2l  1 uv  +  *202v2  +  &330"3  +  &321  +  bH2uv^  +  fe303 v3 


Again,  the  transformations  are  reformulated  in  matrix  form.  Let 
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such  that, 


and 


X  =  WA,  Y  =  WB 


e  =  minA  8[(X- WA)T(X- WA)  +  (Y -  WB)7 (Y  -  WB)] 


(14) 


The  solution  to  Eq.  (14)  is 


A  =  W~]X,  B  =  W  'Y 


t  y  — 1  T 

where  W~  is  the  pseudo-inverse,  computed  by  ( W  W)  W  . 

7.5.2  Thin-Plate  Spline  Interpolation 

Interpolation  differs  from  approximation  in  the  way  that  interpolation  problem  tries  to  construct  a  continu¬ 
ous  function/:  9trt  9tm  so  that f(pi)  =  q where  p\  and  qx  are  finite  sets  of  distinct  point  of  91*  and  9tm 

respectively;  while  approximation  can  only  find  function  g  :  91"  9tm  such  that  #(/?,)  ~  q- .  In  order  to 
achieve  well-posed  (existence,  uniqueness,  and  stability)  solution  to  the  interpolation  problem,  it  proves 
quite  natural  to  require  the  minimization  of  some  appropriate  expression,  such  as  the  quadratic  seminorm1 

|/|2  [Meinguet  84].  |/|2  can  be  physically  interpreted  as  the  bending  energy  of  a  thin  plate  of  infinite  extent. 
Therefore,  interpolants  of  minimum  seminorm  can  be  appropriately  termed  thin  plate  splines  (TPS). 


1.  A  real-valued  functional  \x\ ,  defined  on  a  vector  space  X,  is  called  a  seminorm  if:  a)  \x\  >0  for  all  jc  €  X ; 
b)  \ax\  =  |oc|  •  \x\  ;c)  |*  +  y|<M  +  |y| 
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Bookstein  [Bookstein  78][Bookstein  89]  is  the  first  researcher  to  use  TPS  to  model  the  deformations 
between  two  sets  of  landmark  points  (also  called  control  points),  where  TPS  warps  one  picture  into  the  space 
of  another  by  mapping  landmarks  exactly  onto  landmarks  while  minimizing  a  plausible  bending  energy ,  the 
integral  of  summed  squared  second  derivatives  of  the  mapping  as  Eq.  (15).  In  the  case  that  fix,  y)  is  a  bright¬ 
ness  function,  the  integrand  is  referred  to  as  the  quadratic  variation. 


W2  ■  JJ 


2  ^2 

3/ 

$x2  j 


+  2 


'a2/ 


Kdxdy; 


3/ 


\dxdy 


(15) 


In  the  special  case  of  geometric  correction,  TPS  represents  the  mapping  of  control  points  from  the  corrected 
image  to  their  correspondence  in  the  measured  image  in  the  manner  which  uniquely  minimizes  the  bending 
energy  of  Eq.  (15). 


In  Eq.  (15),  A2/  -  (d2f/dx2)2  +  2 (d2 f/dxdy)  +  (d2  f/dy2)2  is  the  biharmonic  equation ,  and  the  fundamen- 

tal  solution  of  A  /  =  0  is  U(r)  =  r  logr  ,  where  r  is  the  distance  between  two  control  points.  Since  A  /  is 
non-negative,  the  unique  solution  that  minimizes  the  bending  energy  is  a  linear  combination  of  U(r) ,  which 
is  bounded  and  asymptotically  flat.  If  scaling  or  homogeneous  shear  is  allowed  before  applying  the  displace¬ 
ments,  an  affine  transformation  should  be  added  to  the  solution  to  indicate  the  behavior  of  /  at  infinity. 
Therefore,  the  complete  solution  to  Eq.  (15)  is  Eq.  (16): 

m 

/(*>  y)  =  X  wjU(-rj'>  +  a\+  axx  +  v  O6) 

;  =  i 

where  m  is  the  number  of  control  points.  Eq.  (16)  consists  of  two  parts:  the  non-linear  part  which  is  the  sum 
of  functions  U(r) ,  and  the  linear  part  which  is  the  affine  transformation. 


In  our  imaging  system,  we  adopt  TPS  to  correct  geometric  distortion.  We  still  assume  (xit  yt)  to  be  a  set  of 
points  from  the  measured  image,  and  («;,  vz)  from  the  corrected  image.  The  TPS  transformation  function  T : 

R 2  -» R 2  is  sought  such  that 

[■*,•>]  =  [t^K,  V,)  Ty(ujt  V,.)]  (17) 


and  T  minimizes  the  bending  energy  of  Eq.  (15).  Substituting  Eq.  (16)  into  Eq.  (17),  we  get 


m  ~  1 


X  PiU{-rij)  +  a\  +  axui +  ayVj 

7  =  0 

m-  1 


X  qiU(rij)  +  bl+bxui  +  byvi 
J  =  0 


7 


where  r y  =  | (w,.,  v-)  -  («;-,  v;)| .  Again,  we  rewrite  the  problem  in  matrix  form  as  Eq.  (18), 
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which  we  denote  as  Eq.  (19), 


[^/nxm  ^mx3J 


r  m  x  1  ^ m  x  1 
^3x1  ^3x1 


[^fflxl  ^  m  x  l] 


(19) 


To  solve  the  coefficient  matrices  P,  Q ,  A ,  and  P,  matrix  should  be  inverted.  To  make  the  matrix  size 

compatible  for  inverting,  Eq.  (19)  is  reconstructed  into  Eq.  (20),  where  O  is  the  zero  matrix  with  appropriate 
dimensions.  The  coefficients  can  then  be  solved  by  Eq.  (21). 
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P  Q 

X  Y 
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(21) 


Similar  to  the  radiometric  correction:  with  the  control  points  chosen,  we  can  construct  matrices  K ,  L,  X,  and 
Z  and  apply  them  to  Eq.  (21)  to  solve  for  the  coefficient  matrices  P,  Q,  A,  and  B.  With  the  coefficients  deter¬ 
mined,  and  the  corresponding  control  points,  a  look-up  table  can  be  calculated  with  each  entry  indicating  the 
coordinates  that  a  pixel  point  in  the  corrected  image  should  be  looked  up  into  from  the  distorted  image. 

The  more  control  points  chosen,  the  longer  computation  time  needed.  However,  since  this  kind  of  computa¬ 
tion  only  happens  once,  and  the  mapping  can  be  saved  in  a  lookup  table  afterwards,  the  efficiency  should  not 
be  a  concern. 

7.6  Degradation  Analysis 

[Aghdasi  et  al.  94]  has  given  a  detailed  analysis  of  the  degradations  in  an  X-ray  (film)  image  formation  sys¬ 
tem.  The  complete  model  is  illustrated  in  Fig.  41. 


rii  n2  r^3 


Figure  41.  A  model  of  radiographic  image  formation  from  [Aghdasi  et  al.  94]. 

The  input  (gy)  is  the  number  of  X-ray  quanta  received  at  the  screen  per  pixel,  and  the  output  (gg)  is  the 
observed  optical  density  of  each  pixel.  This  model  consists  of  three  sources  of  blur  (hfocalspot9  hscatter  hpm)9 
three  scaling  factors  (in,  rj^,  r|2)>  and  three  sources  of  noise  (/ij,  n 2,  ^3). 

The  cause  of  hjocaispot  is  that  the  focal  spot  on  the  anode  of  the  X-ray  tube  is  not  a  perfect  point  source, 
instead,  it  has  a  finite  size,  which  generates  the  blur  effect  as  shown  in  Fig.  42.  hscatter  results  from  the  X-ray 
scatter  in  the  object,  and  hjnm  from  the  film. 
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Figure  42.  Blur  caused  by  finite-size  X-ray  source. 

Of  the  three  scaling  factors,  m  is  the  mean  screen  amplification  ratio,  rji  and  r(2  are  the  absorption  efficien¬ 
cies  of  the  screen  and  the  film  respectively. 

Noises  mainly  come  from  the  film  grain  noise  (n3)  and  the  quantum  noise  due  to  the  discrete  nature  of  the  X- 
ray  photons.  Quantum  noise  can  be  further  divided  into  the  correlated  components  of  the  input  quantum 
noise  («j),  and  the  uncorrelated  components  of  the  input  quantum  noise  [n^). 

Our  imaging  system  is  different  from  Aghdasi’s  since  we  use  CCD  detectors  instead  of  film.  This  has 
avoided  all  the  degradations  caused  by  film,  but  at  the  same  time  brought  new  types  of  degradations.  For 
example,  the  butted  fiber  optic  tapers  introduce  missing  data  along  the  butting  edge.  Therefore,  a  modified 
degradation  model  need  to  be  formulated. 

Blur.  In  our  CCD  imaging  system,  the  blur  caused  by  hjnm  is  avoided.  We  also  assume  that  the  X-ray  scatter 
can  be  largely  absorbed  by  a  vibrating  grid.  Therefore,  the  major  blur  source  comes  from  the  finite-sized  X- 
ray  source.  A  simple  computation  shows  the  effect  of  blur  from  this  source: 

In  a  digital  mammographic  system,  it  is  reasonable  to  assume  that  the  point  source  is  0.6mm  width,  the  aver¬ 
age  distance  from  an  object  to  the  CCD  detector  is  2cm,  and  the  distance  from  X-ray  source  to  the  object  is 
64cm,  then  the  captured  image  of  the  focal  spot  will  be  18[im.  That  is,  the  blur  is  about  one  pixel. 


0.6  mm 

-H  H— 


X-ray  source 


a  point  in  the 
original  image 


captured  image 


64  cm 


2  cm 


Figure  43.  An  example  of  blur  effect  by  finite  size  focal  spot. 

Noise.  Noise  mainly  comes  from  three  sources  in  our  system:  1)  dark  current  noise  which  can  be  reduced  by 
cooling;  2)  quantum  noise,  a  signal-depend  noise  which  can  be  modeled  by  a  Poisson  process  [Rabbani  88]. 
When  the  readout  rate  is  very  high  and  the  amount  of  the  readout  data  is  very  large,  this  kind  of  noise  can  be 
neglected;  3)  fixed  noise,  due  to  the  electrical  system  noise,  quantization  noise,  etc.,  which  is  the  major  noise 
source  in  our  system. 
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CCD  Detector  Induced  Errors.  If  a  single  CCD  cell  is  defective,  it  can  cause  data  from  the  rest  of  the  row 
missing  due  to  CCDs  serial  readout  mechanism.  Also,  butting  technique  can  cause  data  at  the  boundary  of 
adjacent  CCDs  missing.  Both  of  these  sources  introduce  missing  column(s)/row(s)  in  images.  Defects  in  the 
metallization  or  fibers  can  introduce  single  point  defects  as  well. 

7.7  Regularization  Theory 

Missing  data  estimation  is  basically  an  image  restoration  problem.  Image  restoration  is  well  known  to  be  ill- 
posed  [Andrews  and  Hunt  77][Katsaggelos  91],  since  a  little  change  in  the  input  data  can  dramatically  affect 
the  solution.  Ill-posedness  comes  from  the  ill-conditioning  of  the  matrix  constructed  from  the  system  PSF. 
The  most  popular  solution  to  ill-posedness  is  regularization.  Although  many  different  kinds  of  image  resto¬ 
ration  algorithms  have  been  proposed  so  far,  they  all  share  a  common  structure:  the  regularization  theory. 

Generally  speaking,  any  regularization  method  tries  to  analyze  a  related  well-posed  problem  whose  solution 
approximates  the  original  ill-posed  problem.  The  well-posedness  is  achieved  by  implementing  one  or  more 
of  the  following  basic  ideas  [Nashed  76]: 

•  change  of  the  concept  of  a  solution; 

•  restriction  of  the  data; 

•  change  of  the  space  and/or  topologies; 

•  modification  of  the  operator  itself; 

•  the  concept  of  regularization  operators;  and 

•  well-posed  stochastic  extensions  of  ill-posed  problems. 

According  to  [Tikhonov  and  Arsenin  77],  a  regularization  method  consists  of  finding  regularizing  operators 
that  operate  on  the  data,  and  determining  the  regularization  parameter  from  supplementary  information  per¬ 
taining  to  the  problem.  The  regularization  operator  depends  continuously  on  the  data  and  results  in  the  true 
solution  when  the  regularization  parameter  goes  to  zero,  or  equivalently  when  the  noise  goes  to  zero. 

In  the  traditional  image  restoration  approaches,  image  formation  is  modeled  as  Eq.  (22), 

Hf  +  n  =  g  (22) 

where/,  w,  and  g  are  N  x  1  vectors,  the  lexicographic  representation  of  the  original  image,  noise,  and  the 

measured  image  respectively;  and  H  is  the  N2  x  N2  Toeplitz  matrix,  derived  from  the  system  PSF.  The  regu¬ 
larization  method  constructs  the  solution  as  Eq.  (23), 

minf[u(f ,  g)  +  pv(/)]  (23) 

where  /  is  sought  to  minimize  u(f,  g)  +  pv(/).  The  first  term,  u(f \  g ),  describes  how  the  real  image  data  is 
related  to  the  measured  data.  In  other  words,  this  term  models  the  characteristic  of  the  imaging  system.  The 
second  term  is  the  regularization  term  with  the  regularization  operator  v  operating  on  the  original  image  /, 
and  the  regularization  parameter  p  used  to  tune  up  the  weight  of  the  regularization  term.  By  adding  the  regu¬ 
larization  term,  the  original  ill-posed  problem  turns  into  a  well-posed  one,  that  is,  the  insertion  of  the  regu¬ 
larization  operator  puts  some  constraints  on  what/might  be,  which  makes  the  solution  more  stable. 

In  our  problem,  besides  noise  and  blur,  another  corruption,  missing  data,  is  inserted  to  the  measured  image  g. 
The  consistency  method  using  separable  deblurring  constructs  a  well-posed  problem  by  putting  several 
restrictions  on  the  image  data  and  the  blur  kernel.  The  MAP  method  is  actually  a  Bayesian  interpretation  of 
regularization  problems,  which  is  solved  by  a  global  optimization  technique,  called  mean  field  annealing 
(MFA). 
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7.8  Degradation  and  Distortion  Characterization 

Related  parameters  include  degrees  of  radiometric  and  geometric  distortions,  variances  of  point  spread  func¬ 
tion  and  noise  distribution,  and  the  amount  of  missing  columns. 

7.8.1  Degree  of  Radiometric  Distortion 

For  vignetting-type  radiometric  distortion,  the  brightness  degradation  can  be  modeled  by  a  series  of  concen¬ 
tric  ellipses.  We  define  the  degree  of  radiometric  distortion  (r|r)  by  the  eccentricity  of  the  most  inner  ellipse 
(all  the  concentric  ellipses  have  the  same  eccentricity).  T[r  is  thus  equal  to  ratio  of  the  length  of  the  major  axis 
(a)  to  the  length  of  the  minor  axis  ( b ),  as  rir  =  a/b . 


Figure  44.  A  real  image  used  to  compute  the  degree  of  both  radiometric  distortion  and  geometric  distortion. 


We  use  the  image  in  Fig.  44  to  compute  x\r  of  both  the  left  sensor  image  and  the  right  sensor  image  as  Eq. 
(24), 


=  m  = 117 

T\r(right)  =  ^  =  1.47 


(24) 


7.8.2  Degree  of  Geometric  Distortion 


The  degree  of  geometric  distortion  (v\g,  mainly  pincushion  distortion)  is  measured  at  the  smaller  end  of  the 
fiber-optic  taper.  The  r\g  in  the  long  direction  is  defined  as  Eq.  (25)  by  the  manufacturer,  where  Lc,  Lm,  Ls} 
and  Ls2  are  illustrated  in  Fig.  45. 


Kv  2 


(25) 
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Figure  45.  Parameter  illustration  in  the  definition  of  pincushion  distortion  degree. 


By  obtaining  the  corresponding  Lsh  Ls2,  Lm>  and  Lc  from  both  the  left  and  right  parts  of  the  original  image  in 
Fig.  44,  we  calculate  the  degree  of  pincushion  distortion  as  Eq.  (26). 


\(right)  = 


1075  -  1043 
1078  +  1175 
1086-1049 
1171  +  1180 


x  100  =  1.42% 


x  100  =  1.57% 


(26) 


7.8.3  Point  Spread  Function 

The  point  spread  function  (PSF)  can  be  estimated  by  the  derivative  of  the  edge  spread  function  (ESF).  Fig. 
46  shows  the  ESFs  and  the  corresponding  PSFs  along  both  horizontal  and  vertical  directions  at  different 
locations  (from  both  left  and  right  parts  of  the  image,  select  samples  at  top-left,  top-middle,  top-right,  mid¬ 
dle-left,  middle-middle,  middle-right,  bottom-left,  bottom-middle,  and  bottom-right  respectively)  of  the 
image  in  Fig.  44,  where  the  derivatives  are  taken  numerically. 
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Figure  46.  ESFs  (left  column)  and  PSFs  (right  column)  along  the  horizontal  direction  (top)  and  the  vertical  direction 
(bottom). 


July  15, 1999 


51 


We  observe  that  the  PSFs  are  all  behave  like  a  Gaussian,  and  the  ESFs  are  like  an  error  function.  We  there¬ 
fore  fit  the  sample  data  to  the  error  function  instead  of  the  Gaussian. 

The  best  fitting  standard  deviations  (a)  for  each  of  the  error  functions  in  Fig.  46  are  listed  in  Table  4,  and  fur¬ 
ther  illustrated  in  Fig.  47.  We  observe  that  the  values  of  a  are  very  close  at  different  positions  and  different 
edge  orientations.  In  other  words,  the  shapes  of  different  ellipses  in  Fig.  47  are  very  similar,  which  means  the 
PSF  is  roughly  homogeneous.  Therefore,  we  substitute  G  =  1.85  to  the  2-D  Gaussian  function  (Eq.  (27)),  and 
sample  it  to  create  a  5  x  5  blur  kernel  (h' )  as  Eq.  (28).  Kernel  ti  is  further  normalized  to  h  as  Eq.  (29). 


TABLE  4.  The  best  fitting  standard  deviation  for  ESE 


left 

part  of  the  image 

right  part  of  the  image 

along  the  hori¬ 
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tion 

left 

middle 

left 
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Figure  47.  Illustration  of  the  fitting  a  (aA,  oy)  value  from  Table  4  by  ellipses.  The  length  of  the  axes  of  each  ellipse  is 
taken  from  a  along  the  horizontal  direction,  and  a  along  the  vertical  direction  at  the  same  location  of  the  image.  For 
example,  the  axes  of  the  top-left  ellipse  is  (1.82,  1.83)  taken  from  cells  (1,  1)  and  (3, 1)  of  the  table. 
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0.0506 
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0.0677 
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0.0506 
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7.8.4  Noise  Characterization 


Noise  is  characterized  by  taking  a  homogeneous  segment  from  the  captured  image,  and  observing  its  histo¬ 
gram  distribution.  Fig.  48  shows  a  flat  frame  captured  from  our  imaging  system.  We  select  nine  homoge¬ 
neous  segments  at  different  locations  (top-left,  top-middle,  top-right,  middle-left,  middle-middle,  middle- 
right,  bottom-left,  bottom-middle,  and  bottom-right)  from  both  the  left  and  right  parts  of  the  flat  frame  and 
draw  their  histograms  as  Fig.  49.  From  Fig.  49,  we  observe  that  all  the  histograms  behave  like  Gaussian. 
Again,  we  use  interopt  to  fit  the  histograms  to  Gaussian  and  obtain  the  average  best  fitting  standard  deviation 
which  is  3.9. 


Figure  48.  The  image  of  a  flat  frame. 


140  160  180  200  140  160  180  200 


200  220  240  140  160  180  200  140  160  180  200 


Figure  49.  Noise  distribution  of  a  homogeneous  segment  from  nine  different  locations  of  both  the  left  (left  column) 
and  right  (right  column)  parts  of  the  image. 


This  method  gives  us  a  hypothetical  idea  about  how  noise  contributes  to  the  formation  of  the  image. 
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7.8.5  Missing  Data  Identification 

We  measured  that  the  thinnest  line  our  CCD  sensor  can  detect  is  about  50.8pm  (0.002in),  that  is,  the  width 
of  a  pixel  is  50.8pm.  In  order  to  identify  how  many  columns  are  missing  at  the  butting  edge,  we  take  images 
of  a  Precision  Ronchi  ruling  across  the  butting  edge  and  not  across  the  butting  edge. 

Ronchi  rulings  are  evenly  spaced  lines  running  parallel  to  each  other  with  line  width  equal  to  space  width. 
The  ruling  we  use  has  250  lines  per  inch.  Therefore,  the  line  width  of  the  ruling  can  be  calculated  by  Eq. 
(30).  Since  the  width  of  each  line  in  the  ruling  matches  the  width  of  the  thinnest  line  the  sensor  can  detect, 
the  number  of  missing  columns  can  be  identified  by  counting  the  difference  of  number  of  lines  from  ruling 
image  across  the  butting  edge  and  ruling  image  fully  contained  in  one  sensor,  as  shown  in  Fig.  50. 

„  *tw — -  x  25.4cm  =  0.0508cm  =  50.8[im  (30) 

250lines  X  2  r 


Figure  50.  Images  of  Ronchi  rulings  across  the  butting  edge  and  not  across  the  butting  edge. 

We  count  number  of  lines  across  different  positions  of  the  ruling,  with  the  results  that  the  left  image  detects 
237  lines,  and  the  right  image  detects  238  lines.  Therefore,  we  verify  that  the  butted  CCD  sensor  has  one 
column  missing  at  the  butting  edge. 

7.9  The  Consistency  Method  with  Separable  Deblurring 

By  putting  some  restrictions  on  input  data,  and  the  blur  operator,  we  are  able  to  transform  the  missing  data 
estimation  problem  into  a  well-posed  one.  Several  assumptions  are  made  to  simplify  the  problem.  Possible 
relaxation  of  these  assumptions  are  further  discussed  in  Appendix  7.9.5. 

•  the  blur  kernel  is  separable; 

•  the  blur  kernel  is  exactly  known; 

•  noise  is  not  considered  at  this  stage  of  study; 

•  only  one  column  of  data  is  missed  in  the  original  image;  and 

•  the  original  image  is  of  integer  type. 

The  image  formation  is  modeled  as  Eq.  (31), 

f^>h-  g-n  (31) 

where  g  is  the  measured  image  with  one  column  of  data  missing, /is  the  original  image,  and  n,  the  noise,  is 
regarded  as  small  perturbation  of  g  and  is  neglected  at  this  stage  of  study.  The  convolution  kernel  h  is  dis¬ 
crete  and  of  finite  support.  Both/and  g  are  M  x  N  matrices. 
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The  consistency  method  is  proposed  based  on  the  assumption  that  the  blur  kernel  is  separable,  such  as  Gaus¬ 
sian.  The  separability  property  of  Gaussian  leads  to  the  separable  deblurring  where  two  NxN  sparse  matri¬ 
ces  are  generated  based  on  the  separated  blur  kernels.  We  demonstrate  that  this  transformed  problem  is  well- 
conditioned. 

The  consistency  approach  concerns  three  issues:  separable  deblurring,  Householder  triangularization  based 
QR  factorization  (HHQR),  and  consistency  in  missing  data  estimation  [Qi  et  al.  98]. 


7.9.1  Separable  Deblurring 


Since  the  convolution  kernel  h  is  discrete  and  of  finite  support,  it  can  be  written  as  a  matrix  like  Eq.  (32).  If  h 
is  also  separable,  then  the  convolution  can  be  performed  separately  as  Eq.  (33),  where  hy  and  hx  are  the  sep¬ 
arated  vertical  and  horizontal  components  of  h.  •  denotes  the  matrix  multiplication,  and  <8>  is  the  convolu¬ 
tion  operator. 


h  = 


h  11  •••  h\n 


\Ki  -  K 

f®h  =  /  <8>  (hy  •  hx)  =  (/  ®  hy)  <8>  hx  =  (/  <8>  hx)  <8>  hy 


(32) 


(33) 


Assume/is  anMxiV  image  and  h  is  an  m  x  n  Gaussian  kernel,  the  convolution  of  image/ with  the  horizontal 
kernel  component  hx  can  be  equally  achieved  by  a  matrix  multiplication  between /and  aniVxiV( 


n-  1 


'rt-1 

2 


)-band  matrix  Dx  (Eq.  (34)),  generated  from  hx\  and  the  convolution  with  hy  can  be  similarly 


achieved  by  a  matrix  multiplication  with  an  M  x  M  ( 
from  hy. 


I  m-1  I  [m-  1 

L  2  J’  2 


)-band  matrix  Dy  (Eq.  (35)),  generated 


n-l 

2 


</-;< 


otherwise 


m  -  1 
2 


<i~j< 


otherwise 


(34) 


(35) 


For  example,  if  we  have  an  image/of  dimension  5x4,  and  a  separable  blur  kernel  h  of  dimension  4x3: 


/  = 


, 

/ill  h\ 2  /il3 

hy  1 

,  h  = 

hl\  h22  h23 

— 

hy2 

r  ~ 

1 

h3\  h32  h33 

hy3 

y 

y 

_/i41  hA2  h\3 

hyA 

[hx  1  hx 2  hxi\  > 


the  matrices  Dx  and  Dy  then  look  like: 
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With  Dx  and  Dy,  we  can  rewrite  Eq.  (33)  into  Eq.  (36),  which  gives  us  another  interpretation  on  how  the 
original  image  is  blurred.  Eq.  (37)  is  the  inverse  problem  to  Eq.  (36),  where  if  matrices  Dx  and  Dy  are  both 
well-conditioned,  the  solution  /  should  be  stable.  We  will  analyze  the  conditioning  of  Eq.  (36)  in  detail  in 
Appendix  7.9.4. 


t  t 

g  =  (f( 8>  hx)  <S>  hy  =  Dy  •  (Dx  •  /  ) 


(36) 


/  =  ( Dx~' 


{Dy~l 2 


T  T 


• g )  ) 


(37) 


7.9.2  HHQR 

To  actually  solve  the  system  of  Eq.  (36),  we  use  Householder  triangularization  based  QR  factorization 
(HHQR)  algorithm  [Trefethen  and  Bau  97].  HHQR  has  been  proven  to  be  backward-stable  which  means  if 
the  linear  system  is  well-conditioned,  this  algorithm  can  provide  the  accurate  solution.  Assume  the  linear 
system  is  Ax  =  b ,  HHQR  can  then  be  implemented  in  the  following  three  steps: 

•  using  QR  factorization  to  factor  A  into  Q  •  R  ,  i.e.,  QRx  =  b  where  Q  is  a  unitary  matrix  and  R  is  an  upper 
triangular  matrix  generated  by  HH  triangularization; 

-1  T 

•  since  the  inverse  of  a  unitary  matrix  is  equal  to  its  transpose,  Rx  -  Q  b  =  Q  b  =  y; 

•  since  R  is  an  upper  triangular  matrix,  we  can  solve  x  by  simple  scaler  multiplication  and  subtraction 
without  inverse  operation. 

In  our  problem,  HHQR  algorithm  is  used  twice  to  solve  the  two  linear  systems  derived  from  Eq.  (36): 

(1)  solve  g  =  Dyfy  for  fy ,  and 

(2)  solve  fy  =  Dx*  fT  for/. 

7.9.3  Missing  Data  Estimation  by  Consistency 

Assume  the  original  image  has  only  integer  brightness  values,  and  these  values  are  within  [0,  255],  in  partic¬ 
ular,  the  missing  pixels  should  have  their  original  integer  brightness  values  between  0  and  255,  the  two  linear 
systems  are  solved  as  follows: 

(1)  Solve  g  =  Dy  f  for  fy  using  HHQR.  The  process  is  to  deblur  image  g  along  the  vertical  direction.  We 
show  in  Appendix  7.9.4  that  matrix  Dy  is  well-conditioned  such  that  the  inverse  problem  is  well-posed. 

(2)  For  the  second  linear  system  fTy  -  Dx*  fT ,  HHQR  can  not  be  used  directly,  since  there  is  one  column  of 

missing  pixels  in  /  .  Instead,  we  solve  each  row  of  the  image  separately,  i.e.  solve  fTy(i)  =  Dx  •  fT (i) , 

where  i  represents  the  ith  row  in  that  image.  For  the  example  we  took  at  the  end  of  Appendix  7.9.1,  the  sec¬ 
ond  linear  system  can  be  written  as  Eq.  (38), 
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*/,(U) 

hx 2  hx3  0  0 

/(U) 

fy(i,  2) 

hxl  hx2  hx 3  0 

• 

fd,  2) 

fyd,  3) 

0  hxl  hxl  hx3 

fd,  3) 

/yd,  4) 

0  0  hxl  hx 2_ 

fd,  4) 

(38) 


Assume  the  missing  pixel  is  in  the  3rd  column,  we  then  have  5  unknowns  (jf(/,l),  //,2),  //,3),  //,4)  and 
j/y(/,3))  but  only  4  equations  in  the  linear  system.  To  solve  this  problem,  we  need  to  find  another  condition. 
Since  we  already  know  that//, 3)  should  be  an  integer  between  0  and  255,  an  exhaustive  search  can  be  car¬ 
ried  on  by  assuming//, 3)  is  each  one  of  them,  which  gives  us  4  unknowns  and  4  equations.  Reconstruct  vec¬ 
tors  fy ,/,  and  matrix  Dx  in  Eq.  (38),  so  that  all  the  unknowns  are  at  one  side  of  the  equation,  like  Eq.  (39). 

We  denote  the  matrix  reconstructed  from  Dx  as  Drx ,  and  show  in  Appendix  7.9.4  that  it  is  also  well-condi¬ 
tioned. 


fyd,  1) 

hxl  hx 3  0  0 
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fyd,  2) -hxl  ■  fd,  3) 
-hxl  ■  fd,  3) 
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hxl  hxl  0  0 

0  hxl  -1  hx3 

• 

fd,  2) 

fyd,  3) 

fyd,  4) -hxl  •/(,',  3) 

000  hxl 

JdA). 

Solve  the  new  linear  system  by  HHQR,  and  check  the  solution  to  see  if  it  satisfies  the  consistency  criterion , 
that  is,  the  solutions  (/(/,l),//,2),//,4))  are  all  between  0  and  255  and  are  all  integers.  If  they  do,  then  we 
can  go  on  and  solve  the  next  row;  otherwise,  next  value  of//, 3)  is  attempted.  The  well-conditioned  problem 
and  the  backward-stable  algorithm  assure  that  the  solution  is  unique  and  accurate. 

How  to  choose  the  initial  value  of  the  missing  pixel  plays  an  important  role  in  speeding  up  the  exhaustive 
search.  We  use  the  average  brightness  of  the  missing  pixel’s  left  and  right  neighbors  as  the  first  attempt,  and 
pick  up  the  next  attempt  by  swinging  around  this  value  with  increasing  step  (Fig.  51).  In  practice,  this  saves 
half  of  the  time  compared  to  the  normal  approach  which  always  starts  from  brightness  0. 


4th  2nd  3rd  5  th 

-2-1  +1  +2 

I - 1 - 1 - 1 - 1 - 1 - 1 

0  1st  255 

attempt 

Figure  51.  Choose  initial  value  for  the  missing  pixel  in  exhaustive  searching. 

7.9.4  Conditioning  Analysis 

Conditioning  of  a  mathematical  problem  is  measured  by  the  sensitivity  of  output  to  changes  in  input.  For  a 
well-conditioned  problem,  a  small  change  of  input  does  not  affect  the  output  much;  while  for  an  ill-condi¬ 
tioned  problem,  a  small  change  of  input  can  change  the  output  a  great  deal. 

Condition  number  is  the  measurement  of  the  conditioning  of  a  problem.  Generally,  it  is  defined  as  Eq.  (40). 
The  larger  the  condition  number,  the  more  ill-conditioned  the  problem  is. 

condition  number  =  (40) 

change  in  input 

The  conditioning  of  a  linear  system  Ax  =  b  is  determined  by  the  condition  number  of  matrix  A.  The  relative 
condition  number  K  is  defined  as  Eq.  (41), 

K  =  MllUil  (41) 

where  ||.||  usually  indicates  the  2-norm.  K  is  in  the  range  of  [1,  °°) .  When  K  »  1 ,  the  linear  system  is  ill-con¬ 
ditioned. 
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While  most  of  the  literature  on  image  restoration  concentrates  on  smoothing  methods  to  reduce  the  effect  of 
noise,  few  [Forbes  and  Anh  94][Milinazzo  et  al.  87]  pay  attention  to  quantify  the  amount  of  ill  condition  of 
the  blur  kernel.  Since  the  whole  idea  of  the  consistency  method  is  built  upon  the  well-conditioning  of  the 
two  linear  systems,  analyzing  the  amount  of  ill  condition  of  matrices  becomes  an  important  issue. 

In  the  missing  data  estimation  problem,  the  conditioning  of  the  two  linear  systems  is  measured  by  the  condi¬ 
tion  number  of  matrices  Dx ,  Dy,  and  Drx.  Since  Dx  and  Dy  has  similar  structure,  we  only  analyze  Dy  and 
Drx. 

Condition  number  of  Dy.  We  compute  the  condition  number  of  matrix  Dy  generated  with  different  image 
sizes  and  Gaussian  blur  kernel  sizes.  The  results  are  listed  in  Table  5.  We  can  see  that  for  kernel  size  larger 
than  3,  the  condition  number  does  not  change  much  when  increasing  the  image  size.  Also,  when  increasing 
the  kernel  size,  the  condition  number  always  fluctuates  around  10. 


TABLE  5.  Condition  number  of  matrices  Dy. 


Image  Size  (N 

x  AO 

Condition  Number  of  Different  Kernel  Sizes  (n  x  n) 

3 

5 

7 

9 

11 

13 

15 

48.3742 

9.9852 

9.8348 

8.2826 

8.6213 

8.2726 

100 

4133.6000 

9.9659 

10.9641 

8.9596 

9.3753 

8.4998 

1000 

406100.0000 

9.9996 

10.9996 

8.9795 

9.3997 

8.8548 

2000 

9.9999 

10.9999 

8.9797 

9.3999 

8.8550 

Fig.  52  is  several  plots  illustrating  the  rate  of  growth  of  the  condition  number  with  respect  to  different  image 
sizes  and  kernel  sizes.  Fig.  52  (a)  is  the  growth  curve  for  3  x  3  blur  kernel.  The  curve  increases  steadily  with¬ 
out  converging  to  an  upper  bound,  which  indicates  the  system  is  ill-conditioned. 

Fig.  52  (c)  shows  six  curves  with  respect  to  six  different  kernel  sizes:  5x5, 7x7,  9x9,  11x11,13x13,  and 
15  x  15.  All  the  curves  behave  similar:  a  sharp  increase  at  the  beginning  and  asymptotically  converge  to  a 
constant  (less  than  11)  after  the  image  size  is  larger  than  100  x  100.  Fig.  52  (b)  displays  three  curves  corre¬ 
sponding  to  different  image  sizes.  It  indicates  that  all  the  condition  numbers  fall  into  the  area  with  an  upper 
bound  created  by  the  largest  image  size  (2000  x  2000)  and  a  lower  bound  related  to  the  smallest  image  size 
(15  x  15).  The  upper  bound  is  around  10,  which  demonstrates  the  well-conditioning  of  the  system. 

From  the  above  analysis,  we  claim  that  when  blur  kernel  is  larger  than  3x3,  the  first  linear  system 
g  =  Dy  fy  is  well-conditioned. 


Figure  52.  Condition  number  of  matrix  Dy  with  respect  to  different  image  and  kernel  sizes. 

Condition  number  of  Drx.  The  condition  number  of  Drx  indicates  the  conditioning  of  the  second  linear 
system  (Eq.  (39)).  Drx  is  transformed  from  Dx ,  with  only  one  column  different.  If  the  missing  column  in  the 
original  image  is  column y,  then  Drx  is  constructed  by  setting  the  yth  column  of  Dx  to  zero  except  element  (/, 
j)  with  a  value  -1,  as  indicated  in  Eq.  (42),  where  j  =  3. 
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Figure  53.  Condition  number  of  matrix  Drx  with  respect  to  different  image  and  kernel  sizes. 

Fig.  53  is  three  plots  similar  to  Fig.  52.  For  kernel  size  larger  than  3x3,  the  condition  number  still  converges 
asymptotically  to  a  constant  after  an  apparent  increase  at  the  beginning.  Although  the  constant  is  larger  (less 
than  450)  than  that  from  Fig.  52,  it  still  indicates  well-conditioning.  Fig.  53  (b,  c)  also  tells  us  that  the  condi¬ 
tion  number  is  more  affected  by  different  kernel  size  than  different  image  size,  which  can  be  seen  from  the 
close  distance  between  the  three  curves  in  Fig.  53  (b)  and  the  comparatively  larger  distance  between  curves 
in  Fig.  53  (c). 


7.9.5  Assumption  Relaxation  Analysis 

At  the  beginning  of  Appendix  7.9,  we  made  several  assumptions  to  simplify  the  restoration  problem.  This 
section  discusses  the  possible  relaxations  to  these  assumptions.  Except  that  the  blur  kernel  still  needs  to  be 
separable,  we  show  that  other  assumptions  can  all  be  relaxed  to  a  certain  degree,  though  by  sacrificing  cer¬ 
tain  amount  of  accuracy  of  the  solution.  We  first  analyze  sensitivity  of  solution  to  perturbations  in  blur  ker¬ 
nel  ( h ),  and  perturbations  in  measured  image  (g)  caused  by  the  insertion  of  noise  ( n ).  We  also  evaluate  the 
conditioning  of  problem  when  there  are  more  than  one  column  of  missing  data.  Finally,  we  relax  the 
assumption  of  integer-typed  original  image,  and  design  a  neighbor  least  square  error  criterion  to  select  the 
optimal  solution. 

Sensitivity  of  Solution  to  Perturbations  in  Blur  Kernel.  According  to  [Trefethen  and  Bau  97],  in  a  linear 
system  Ax  =  b ,  sensitivity  of  solution  x  to  perturbations  in  b  can  be  measured  by  Eq.  (43), 

I* -ill  _  K(A}\M  (W 

irar  ”  ii  w  {) 


where  tj  =  .  Sensitivity  of  x  to  perturbations  in  A  can  be  estimated  by  Eq.  (44). 

\\Ax\\ 


I* -ill  ^  J8A] 

INI  ~K{A)\\a\\ 


(44) 


In  our  system,  perturbations  in  blur  kernel  ( h )  are  reflected  in  blur  matrices  Dy  and  Drx  as  shown  in  Eq.  (45) 
and  Eq.  (46). 


(Z>y  +  5Dy)/y  =  * 

(45) 

(. Drx  +  8Drx)f  =  f] 

(46) 
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Based  on  Eq.  (43)  and  Eq.  (44),  sensitivity  of/with  respect  to  perturbations  in  Dx  can  be  expressed  by  Eq. 
(47),  and  sensitivity  of/to  perturbations  in  Dy  by  Eq.  (48), 


I1/-/I1 

ii/i 


<  K(Drx)  ■ 


|8ggfl 

\\Drx\\ 


(47) 


\\f-f\\  _K(Drx)  |S/J  lODrx)K(Dy)  |8Dy| 
ll/ll  *1*  '  ||/J  "  r,,  '  \\Dy\\ 


where  7]  x  =  ,  a  number  greater  than  1,  but  very  close  to  1.  If  Drx  and  Dy  have  same  degree  of  per- 

|| UYX  *  /  jj 

turbation,  then  Eq.  (47)  derives  a  smaller  upper  bound  than  Eq.  (48).  Therefore,  sensitivity  of/  to  perturba¬ 
tions  in  blur  kernel  h  is  measured  more  accurate  by  Eq.  (47). 


If  K(Drx) «  10c ,  and  \\&Drx\\/\\Drx\\ «  10  d ,  when  c-d ,  sensitivity  of  solution  is  close  to  1  (Eq.  (49)),  which 
means  the  problem  is  ill-conditioned,  since  we  expect  no  correct  digits  from  the  output/. 


II/- 71  _  m 

i/i  "  ii/i 


<i 


(49) 


Sensitivity  of  Solution  to  Perturbations  in  Measured  Image.  Perturbations  in  measured  image  g  can  be 
interpreted  as  the  effect  of  inserted  noise.  The  two  perturbed  linear  systems  can  be  written  as  Eq.  (50)  and 
Eq.  (51). 


Dy  ■  fy  =  g  +  8g  (50) 

Drx'f  =  fy  +  §fy  (51) 


Sensitivity  of/v  to  perturbations  in  measured  image  g  is  defined  by  Eq.  (52)  based  on  Eq.  (43), 

\\fydi  -  IlSgll 
\\fy\\  h,,  llgll 


(52) 


where  r|v  =  ,  a  number  greater  than  1,  but  very  close  to  1.  Sensitivity  of/to  perturbations  in  mea- 

'  J  y  | 


sured  image  g  can  be  derived  by  Eq.  (53). 


11/ -711  _  KjDrx )  |8/y|  _  K(Drx)K(Dy)  ||8g|| 


ll/ll 


WfJi 


M 


(53) 


Compare  Eq.  (53)  with  Eq.  (47),  if  Drx  and  g  have  same  degree  of  perturbation,  then  solution /can  be  less 
sensitive  to  perturbations  in  blur  kernel  than  to  perturbations  in  measured  image. 

Conditioning  Analysis  for  More  Than  One  Missing  Column.  If  there  are  more  than  one  column  of  miss¬ 
ing  data  in  the  measured  image,  condition  number  of  the  vertical  blur  matrix  Dy  stays  the  same,  but  condi¬ 
tion  number  of  the  reconstructed  horizontal  blur  matrix  Drx  is  changed.  Eq.  (54)  shows  a  template  of  Drx 
transformed  from  Dx  when  two  columns  of  data  (the  2nd  and  the  3rd  columns)  are  missed.  The  missed  ith 
andjth  columns  of  Dx  are  set  to  zeros  except  at  element  (/,  i)  and  (jj)  with  a  value  -1,  similar  to  the  forma¬ 
tion  in  Eq.  (42). 


hxl  hx 3  0  0 

hx  1  hxl  hx3  0 
0  hx\  hxl  hx 3 
0  0  hxl  hxl_ 


hxl  000 
hxl  -100 
0  0-10 
0  0  0  hxl_ 


(54) 
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With  more  than  one  column  missing,  two  problems  need  to  be  considered:  how  many  missing  columns  can 
be  in  the  image,  and  how  far  away  these  missing  columns  should  be  apart  for  stable  inverse.  We  still  use  con¬ 
dition  number  to  measure  the  conditioning  of  the  problem. 

Fig.  54  shows  the  condition  number  of  Drx  with  blur  kernel  ( h )  at  dimension  5x5.  The  three  plots  are  gen¬ 
erated  with  different  numbers  of  missing  columns  (2,  3,  and  4).  We  can  see  that  no  matter  how  many  col¬ 
umns  are  missed  in  the  measured  image,  and  how  far  away  the  two  nearest  missing  columns  are  apart,  all  of 
the  plots  behave  similarly  -  after  an  initial  perturbation,  all  plots  converge  to  a  constant  condition  number, 
which  is  around  120. 


Figure  54.  Condition  number  of  Drx  with  different  numbers  of  missing  columns  at  kernel  size  5x5. 

Fig.  55  shows  the  condition  number  of  Drx  with  blur  kernel  ( h )  at  dimension  7x7.  The  three  plots  also  have 
similar  characteristics  as  those  from  Fig.  54.  The  difference  lies  in  that  after  an  initial  perturbation,  the  three 
plots  in  Fig.  55  converge  to  a  larger  constant  condition  number,  which  is  around  260. 


Figure  55.  Condition  number  of  Drx  with  different  numbers  of  missing  columns  at  kernel  size  7x7. 

Fig.  56  compares  the  condition  number  of  Drx  with  different  kernel  size  when  two  columns  of  data  are 
missed.  Both  of  the  two  curves  increase  steadily  with  respect  to  the  kernel  size. 


Figure  56.  Condition  number  comparison  of  Drx ,  missing  2  columns  of  data,  with  different  kernel  size:  5  x  5, 7  x  7, 
9x9,  and  11x11. 
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These  results  substantiate  our  claim  that  when  more  than  one  column  of  data  are  missed,  it  is  not  the  image 
size,  or  the  distance  between  two  nearest  missing  columns,  or  the  number  of  missing  columns  that  affect  the 
condition  number,  it  is  the  kernel  size  that  plays  the  most  important  role.  When  the  kernel  size  is  less  than  10 
x  10,  the  problem  is  still  well-conditioned.  The  larger  the  kernel  size,  the  worse  the  problem  is  conditioned. 

Neighbor  Least  Square  Error  Consistency  Criterion.  The  last  assumption  we  made  is  that  the  original 
image  needs  to  be  of  integer  type.  We  relax  this  assumption  to  float-type  image  and  design  a  new  consis¬ 
tency  criterion  to  be  used  when  selecting  the  optimal  solution  at  the  last  step  of  the  consistency  method.  We 
call  it  neighbor  least  square  error  consistency  criterion.  It  works  in  the  following  steps: 

1)  Estimate  the  missing  column  in  the  measured  image  using  average  value  of  its  left  and  right  neighbors; 

2)  Apply  HHQR  to  solve  the  first  linear  system  g  =  Dy  f  for  fy\ 

3)  Reconstruct  horizontal  blur  matrix  Dx  to  Drx ,  and  solve  the  second  linear  system  of  Eq.  (39)  one  row  at  a 
time. 

4)  For  each  possible  solution  of  /(/)  (the  zth  row  of  the  restored  image),  compute  the  neighbor  least  square 
error  by  Eq.  (55),  where  d  is  the  column  position  of  the  missing  data,  and  5  is  the  number  of  neighbors  that 
is  involved  in  this  least  square  error  computation.  It  is  important  to  choose  close  neighbors  instead  of  the 
entire  row  for  the  computation. 

</  +  8 

\\f(i)®hx-fy(i)\2  =  X  mi)®hx)j-fy(i,j)f  (55) 

j-  d~  5 

This  criterion  relaxed  the  integer-type  assumption,  unfortunately,  the  search  in  the  fourth  step  has  to  be 
exhaustive.  That  is,  we  need  to  assume  the  value  of  the  missing  pixel  to  be  from  0  to  255.  Based  on  the  256 
solution,  we  select  the  one  that  minimizes  Eq.  (55). 


7.10  MAP  Using  MFA 


This  is  another  approach  we  proposed  to  solve  the  problem  of  missing  data  estimation  with  denoising  and 
deblurring.  This  approach  uses  the  optimization  technique,  MFA,  to  solve  an  MAP  estimate.  With  the  image 
formulated  as  Eq.  (24),  the  maximum  a-posteriori  probability  (MAP)  approach  tries  to  find  an  estimate  of 
image/that  maximizes  the  a-posteriori  probability  p(f\g)  as  Eq.  (56). 


/  =  axgmaxj-p(f\g) 


According  to  Bayes’  rule,  p(f\g)  can  be  written  as  Eq.  (57), 


P(f\g)  = 


p(g\f)P(f) 
P(8 ) 


(56) 


(57) 


where  P(f)  is  the  a-priori  probability  of  the  unknown  image/,  P(g)  is  the  probability  of  g  which  is  a  constant 
when  g  is  given,  and  p{g\f)  is  the  conditional  probability  density  function  (pdf)  of  the  observation  g.  In 
[Szeliski  89],  P(f)  is  called  the  prior  model ,  and  p(g\f) ,  the  sensor  model ,  which  is  a  description  of  the 
noisy  or  stochastic  processes  that  relate  the  original  unknown  image/to  the  measured  image  g. 


Based  on  Bayes’  rule,  the  MAP  approach  can  then  be  derived  as  Eq.  (58),  where  we  neglect  the  constant 

p(gy 

f  =  Mgmaxfp(f\g)  =  zrgmaxf(p(g\f)P{f))  (58) 

The  MAP  approach  may  be  regarded  as  a  Bayes  interpretation  to  the  regularization  theory  [Demoment 
89]  [Li  95].  It  eventually  formulates  an  image  restoration  problem  into  finding  the  optimal  solution  of  an 


July  15, 1999 


62 


objective  function  which  is  very  similar  to  the  traditional  regularization  method,  but  derived  from  a  different 
point  of  view. 

In  the  next  four  subsections,  we  will  discuss  our  design  of  the  sensor  model  and  the  prior  model.  The  mean 
field  annealing  global  optimization  technique  is  adopted  here  to  obtain  the  restored  image  /  which  maxi¬ 
mizes  the  a-posteriori  probability. 


7.10.1  Sensor  Model 


If  we  define  the  5-function  5(0  as  Eq.  (59),  then  the  sensor  model  can  be  expressed  as  Eq.  (60), 


6(0  = 


1 

u 


ie  Q.-Q.J 
is  Qd 


(59) 


n  =  {«/(«/  =  8(0 -[(/ ®  ^  -  g,.] } 


(60) 


where  i  denotes  individual  pixel,  d  denotes  locations  of  missing  data,  £1  is  the  set  of  pixels  of  the  entire 
image,  and  £ld  is  the  set  of  missing  pixels. 


2 

Assume  noise  is  independent  Gaussian  noise  with  zero  mean,  i.e.,  nt~  N(0,  a  ) ,  then  the  likelihood  density 
function  can  be  written  as  Eq.  (61),  a  product  over  the  pixels  of  the  image. 


P(g\f)  = 


n  — 

ieh- n,V2Sa 


exp 


2\ 


2a 


n 


J2na 


(61) 


7.10.2  Prior  Model 

The  prior  model  is  a  probability  description  of  the  original  image.  It  plays  the  same  role  as  the  regularization 
term:  by  providing  some  prior  knowledge  of  the  original  image,  so  to  put  some  constraints  in  the  solution, 
the  original  ill-posed  image  restoration  problem  turns  into  a  well-posed  one. 

The  prior  knowledge  of  the  original  image  refers  to  the  a-priori  belief  that  the  state  of  a  pixel  is  entirely 
determined  by  the  states  of  its  neighboring  pixels.  Specifically,  it  is  expected  that  pixels  close  to  each  other 
tend  to  have  the  same  or  similar  brightness  values  [Besag  74][Besag  86].  [Geman  and  Geman  84]  uses 
Markov  Random  Field  (MRF)  to  represent  the  local  property  of  the  image. 

An  MRF  is  a  probabilistic  process  in  which  all  interaction  is  local.  It  is  defined  over  a  discrete  field  where 
the  probability  of  a  particular  variable^  defends  only  on  a  small  number  of  its  neighbors,  which  is  expressed 
as  Eq.  (62). 

=  P(fi\UjJz  «,-})  (62) 

Although  an  MRF  is  the  correct  model  to  represent  the  local  property  in  the  image,  it  is  difficult  to  estimate 
either  the  conditional  Markov  distribution  p(/,|/) ,  or  the  joint  probability  directly  from  the  conditional  dis¬ 
tribution.  Fortunately,  [Hammersley  and  Handscomb  64]  and  later  paper  [Besag  74]  have  shown  the  equiva¬ 
lence  between  Gibbs  distributions  and  MRF.  This  equivalence  allows  the  modeling  of  local  structure  through 
energies  which  describe  the  interactions  of  pixels  within  each  clique  of  the  neighborhood.  The  a-priori  prob¬ 
ability  of  an  image  by  a  Gibbs  distribution  is  defined  as  Eq.  (63), 

/>(/)  =  exp(-t/(/)/D  (63) 

where  Z  =  ^exp  (-U(f)/T)  is  a  normalizing  constant,  called  the  partition  function’,  T  is  the  temperature  of 
/ 

the  model;  and  U(f)  is  the  energy  function,  that  can  be  written  as  Eq.  (64), 
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utf)  =  %v4f)  (64) 

i 

where  Q  is  the  clique  formed  involving  pixel  i.  Vc{  f)  is  called  the  potential ,  which  depends  on  the  local 
configuration  of  clique  Cv 

The  prior  energy  function  is  usually  formulated  based  on  the  smoothness  property  of  the  original  image. 
Since  the  more  probable  configurations  are  those  with  higher  P(f)  (thus  lower  £/(/)),  U(f)  should  measure  the 
extent  to  which  the  smoothness  is  violated  [Li  95].  For  spatially  continuous  MRFs,  the  energy  function  often 
involves  derivatives.  Different  orders  of  derivatives  imply  different  classes  of  smoothness.  Based  on  the  dis¬ 
cussions  in  [Li  95],  we  summarize  the  different  formulations  of  the  energy  function  according  to  different 
kinds  of  image  surfaces  in  Table  6: 


TABLE  6.  Prior  energy  for  different  surfaces. 


Surface  property 

Derivative 

Prior  energy  U(j) 

flat  surface: 

/(■*.  y)  =  «o 

zero  lst-order 
derivative: 

/,  =  0./v  =  0 

jl(.f2x  +  f2y)dxdy 

mem¬ 

brane 

planar  surface: 
f(x,y)  =  a0  +  alx  +  a2y 

zero  2nd-order 
derivative: 

fxx  =  0,f„  =  o 
fxy  =  0 

quadratic  variation: 
ttifIx+lfxy  +  flyWxdy 

plate 

square  Laplacian: 
j\(fxx  +  fyy)2dxdy 

quadratic  surface: 

2  2 
f(x,y)  =  a0  +  a,jc  +a2xy  +  a2y 

zero  3rd- order 
derivative: 

fxxx  =  °>  fyyy  =  0 
fxx,  =  °’  f xyy  =  0 

jJ(/L+3/L+3/L+/L)^ 

\\(f2xxx  +  f)yy)dxdy 

[Blake  and  Zisserman  87]  designed  a  clipped  parabola  (Fig.  57)  to  model  the  energy  of  interaction  between 
neighbors  in  the  weak  string  (the  1-D  case  of  membrane),  where  V  represents  any  operator  which  returns  a 
measure  of  the  local  “edginess”  of  the  image,  such  as  the  square  Laplacian.  The  central  dip  of  this  energy 
function  punishes  the  difference  Vf(/) ,  and  the  plateaus  allow  discontinuity.  Since  both  noise  and  edges  can 
increase  the  brightness  difference  between  neighbor  pixels,  and  the  model  should  only  punish  the  noise  but 
not  the  edges,  the  energy  curve  turns  into  a  constant  after  V  ■(/)  surpasses  a  certain  threshold.  This  model 

is  in  general  not  convex,  therefore,  graduated  nonconvexity  (GNU)  algorithm  is  developed  which  approxi¬ 
mates  the  energy  with  a  piecewise  smooth  function  as  Eq.  (65), 


Figure  57.  Prior  energy  of  the  GNC  algorithm. 
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V(t)  = 


(65) 


X2t 2 

a-c\\t\-r)2/2 


if  (U|  <q) 
if  ( q  <  |/|  <  r) 
if  (Ul  >  r) 


where  2  =  af-,  +  4:1 ,  and  q  =  4r- . 

XJ  Xr 

In  our  system,  we  adopt  the  energy  model  proposed  in  [Bilbro  and  Snyder  88a],  where  the  penalty  function 
of  Fig.  57  is  interpreted  as  a  Gaussian  with  the  upside  down,  shown  in  Fig.  58.  The  prior  energy  function  can 
then  be  written  as  Eq.  (66), 


Figure  58.  An  ideal  curve  to  represent  energy  function  for  the  prior  model  (the  horizontal  axis  is  the  difference  in 
brightness  of  adjacent  pixels). 


«/>  - 


IV*/,) 

2t2 


(66) 


k 

where  p  is  the  parameter  used  to  adjust  how  smooth  the  image  goes,  V  f.  is  some  norm  of  the  k- th  deriva¬ 
tive,  and  depends  upon  the  property  of  the  image.  The  derivative  operation  can  also  be  formulated  into  a  con¬ 
volution  as  (/®r)2,  with  appropriate  convolution  kernel  r.  Therefore,  the  a-priori  probability  can  be 
written  as  Eq.  (67). 


P(f)  =  z-' 


V®r)2  Y) 

2t2  J 


V 


J 


(67) 


7.10.3  Objective  Function 

With  both  p(g\f)  and  P(f)  developed,  we  can  derive  our  objective  function  as  Eq.  (68)  from  the  MAP  esti¬ 
mate  by  taking  the  natural  logarithms  of  Eq.  (58);  and  by  changing  the  sign,  we  convert  the  problem  from 
maximizing  a  probability  to  minimizing  an  objective  function. 

H  =  -\*(p(g\f)P(f))  =  -\n(p(g\f)-\nP(f))  =  Hn  +  Hp  (68) 


where 
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Hn  =  7~2  X  U/®  *)/-*/] 


2a  is  £j-n 


^  P  f  (/®r)fl 

"  -?A-<exp[“J  (70) 

//„  is  called  the  nowe  term  (Eq.  (69))  which  concerns  both  the  measured  image  g  and  the  original  image /.  Hp 
is  the  prior  term  (Eq.  (70)),  which  is  determined  only  by/.  A  global  optimization  technique  called  mean  field 
annealing  is  used  to  find/which  minimizes  the  objective  function  H  globally. 


7.10.4  Optimization  by  Mean  Field  Annealing 

Mean  field  annealing  (MFA)  is  used  to  find  image  /which  minimizes  the  objective  function  H.  MFA  was 
introduced  in  late  80s  by  different  groups  working  independently  [Bilbro  and  Snyder  88a][Bilbro  et  al. 
89][Peterson  and  Soderberg  89].  It  is  derived  based  on  the  principles  of  simulated  annealing  (SA).  SA  is  able 
to  locate  the  global  minimum/maximum  instead  of  being  trapped  in  local  minima/maxima  because  it  allows 
uphill  moves  with  a  certain  probability.  Both  MFA  and  SA  relate  to  the  physical  process  of  annealing,  where 
at  high  temperatures,  the  particles  in  an  object  are  more  likely  to  randomly  change  their  states;  when  temper¬ 
ature  is  gradually  reduced,  particles  remain  comparatively  stable,  and  a  minimum  energy  state  is  achieved  at 
last.  The  difference  between  MFA  and  SA  is  that  S A  stochastically  simulates  this  process  by  random  search, 
while  MFA  analytically  approximates  this  process  with  a  series  of  deterministic  gradient  descents  [Bilbro  et 
al.  92].  Therefore,  MFA  is  much  more  efficient  than  SA,  as  much  as  50  times  faster  [Bilbro  et  al.  89]. 

Since  its  introduction,  MFA  has  found  its  applications  in  restoration  of  locally-homogeneous  images  [Hiriy- 
annaiah  et  al.  89],  range  images  [Bilbro  and  Snyder  89],  and  locally  smooth  images  [Bilbro  and  Snyder  90]; 
in  image  segmentation  [Snyder  et  al.  91];  in  motion  analysis  [Abdelqader  et  al.  92];  in  sensor  fusion  [Bilbro 
and  Snyder  88b];  in  image  resolution  increasing  [Wang  96],  and  in  solving  anisotropic  diffusion  problems 
[Qi  et  al.  97]. 

In  order  to  use  MFA  in  our  problem,  the  derivative  of  H  with  respect  to /in  the  gradient  descent  equation 
(Eq.  (71))  need  to  be  derived.  Eq.  (72)  shows  the  derivation. 

/+1=/-a|y  (71) 

dH  .  dHp 


a/i  a/i  a/,' 


(n  ®  h  v). 


(3  (/Or) 

— - —exp 


(f®ry 

2t2 


where  hrev  and  are  the  reversed  kernel  of  h  and  r.  n  is  defined  in  Eq.  (60)  where  nx  is  zero  for  unmeasured 
elements.  The  reverse  convolution  with  hrev  will  be  discussed  in  detail  in  Appendix  7.10.5. 

MFA  is  implemented  by  replacing  parameter  %  in  the  Gibbs  distribution  with  X  +  7,  where  7  is  called  the  tem¬ 
perature.  Initially,  7  is  very  large  which  results  in  a  convex  function  minimized  by  /  =  g Gradient  descent 
uses  this  point  as  the  starting  point.  By  gradually  reducing  7,  and  making  it  approach  zero,  a  minimum 
energy  state  will  finally  be  achieved  which  is  where  the  optimal  solution  locates. 


7.10.5  Missing  Data  Estimation  by  MFA 

Compare  to  the  consistency  method,  MFA  is  more  flexible,  less  vulnerable  to  large  amount  of  noise  and  the 
inaccuracy  calculation  of  blur  kernel.  Three  issues  need  to  be  discussed  in  order  to  implement  this  approach: 
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(1)  the  choice  of  the  blur  kernel  (h)\  (2)  the  choice  of  the  neighborhood  operation  kernel  (r);  and  (3)  the 
reverse  convolution  with  hrev  in  the  noise  term. 

The  blur  kernel  ( h ).  How  to  compute  the  blur  kernel  in  an  optimal  way  is  beyond  the  scope  of  this  research. 
[Hussain  97]  states  in  detail  about  how  to  estimate  the  point  spread  function  with  subpixel  accuracy.  Here, 
we  use  the  edge  spread  function  (ESF)  to  estimate  the  point  spread  function  (PSF).  We  first  choose  several 
sets  of  sample  pixels  around  a  horizontal  edge  and  a  vertical  edge  respectively  from  the  measured  image. 
Averaging  the  brightness  of  the  sets  of  sample  pixels  to  obtain  a  discrete  data  set  of  ESF  in  the  horizontal 

and  vertical  direction,  denoted  as  j/^-jy  e  [/-5,  ev  +  5]  j,  and  j  e  [eh-§,  eh  +  8]j,  where  e  and  eh 

are  the  locations  of  the  vertical  and  horizontal  edges  from  the  measured  image,  8  is  a  small  distance  away 
from  the  edge  location.  The  data  sets  {/*;}  and  {/J  .}  should  behave  like  Fig.  59  (a),  shape  of  an  error 

function.  The  derivatives  (or  difference  in  discrete  case)  of  {/}j}  and  {f]j}  (denoted  as  {df^j}  and 
{df]  j] )  should  behave  like  Fig.  59  (b),  shape  of  a  Gaussian. 


Pbetooslion 


Figure  59.  The  edge  spread  function  and  the  point  spread  function. 

The  PSF  is  characterized  by  fitting  the  discrete  data  sets  {d/jj}  and  [df]  j)  to  a  Gaussian,  thus  to  obtain 
the  standard  deviations.  Based  on  the  relationship  between  an  error  function  and  a  Gaussian  function 
(expressed  in  Eq.  (73)),  we  can  instead  fit  the  error  function  directly  with  the  sampled  data  sets  {/?  •}  and 

{/V  j}  .  By  doing  so,  we  avoid  the  difference  operation  and  achieve  higher  accuracy  in  parameter  estimation. 


With  the  standard  deviations  solved,  we  can  then  obtain  the  2-D  PSF,  which  is  shown  in  Eq.  (74).  The  con¬ 
tinuous  PSF  is  sampled  to  construct  a  5  x  5  blur  kernel  h.  Here,  we  use  a  software  package  called  interopt 
[Bilbro  and  Snyder  91]  to  fit  the  sample  data  to  an  error  function. 


G{x ,  y) 


2  ™X<*y 


2  2  2  2 
oyx  +  oxy 


The  neighborhood  operation  kernel  (r).  We  choose  quadratic  variation  (QV)  to  model  the  prior  energy 

2  2 

function.  Compared  to  squared  Laplacian,  QV  does  not  go  to  zero  when  -d-  =  — L  (a  saddle  point)  and 

dx2  dy2 
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conclude  that  a  saddle  point  is  a  plane  as  Laplacian  does.  QV  has  a  more  stable  description  about  the  surface, 
it  also  simulates  the  second  derivative  operation,  and  permits  piecewise  linear  solutions. 

The  discrete  form  of  quadratic  variation  can  be  expressed  as  Eq.  (75), 

J)  +  2f2xyV’  J)  +  f). y(f.  N  (75) 

itj 


where  /„  =  /u_,  +  -2fu,  fxy  =  u+i  +  /;,;+!  +/,-w<  and 

/  =  fi_l  j  +  /f+  j  j  -  2fi  j  [Blake  and  Zisserman  87].  All  of  the  three  2nd  derivatives  can  be  constructed 

into  a  convolution  operation  with  certain  kernels  h ^  and  hyy  as  Eq.  (76). 


2  l].  hxy  = 

0  1 

0 

1 

0-1 

1 

’  hyy  = 

-2 

0  0 

-1 

1 

(76) 


The  reverse  convolution  in  noise  term.  In  Eq.  (72),  the  convolution  with  reverse  kernel  ( hrev )  is  different 
from  others  because  g  has  certain  amounts  of  data  missed.  The  convolution  should  be  taken  at  places  where 
a  pixel  has  a  measurement.  Snyder  derived  the  formulation  for  reverse  convolution  in  [Snyder  99],  which  is 
briefly  described  as  follows: 


Take  the  1-D  example,  assume  the  measured  image  has  every  other  column  of  data  missing  as  indicated  in 
Fig.  60,  where ^  is  a  pixel  from  the  original  image,  gt  is  a  pixel  from  the  measured  image,  and  h  is  the  hori¬ 
zontal  blur  kernel  with  a  finite  kernel  size  5.  We  explain  the  derivative  of  the  noise  term  (Eq.  (72))  in  gradi¬ 
ent  descent  in  detail. 
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Tigure  60.  A  1-D  example  with  every  other  sample  missing  in  the  measured  signal. 

We  first  write  out  all  the  terms  involving  a  pixel  (f4)  at  which  a  measurement  (g4)  was  made  in  the  noise  term 
Hn  according  to  Eq.  (69)  as  Eq.  (77), 

£4  =  ((/  ®  hh  - S2)2  +  ((/  ®  *)4  -  grf  +  ((/  ®  h)6  -  g6 )2 

=  (/o^-2  +  A^-1  +  fl^0  +  +/V*2-S2)  (Jl) 

+  {f  2^-2  +  f  3^-1  +  f  Ah 0  +  f  +  f$h2  ~  g4) 

+  if 4^-2  +  f  1  +  f  0  +  fl^\  +  fs^2 


The  derivative  of  Hn  with  respect  to  pixel /4  can  then  be  derived  as  Eq.  (78),  and  further  generalized  as  Eq. 
(79), 
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(78) 


r  =  2((/®/i)2-S2)/i2  +  2((/®/i)4-«4)/!0  +  2((/®/!)6-S6)/*-2 


t/JI  „ 

=  ((/®A)2-*2)*2  +  O-A1  +  ((/®fc)4-*4)Ao  +  0-A_1 
4  (79) 

+  ((/®A)6-*6>*-2 
=  (((/®*)-«)®A„v)4 

where  /ir£,v  =  h2  h{  h0  h.{  h_2,  and  ( f®h-g )  is  computed  at  all  points  where  g;  is  measured  (that  is,  at  g2, 
and  g6). 

We  also  write  out  all  the  terms  involving  an  unmeasured  pixel  (say/3)  in  the  noise  term  Hn  as  Eq.  (80). 

E3  =  ((/  ®  h)2  -  g2)2  +  ((/  ®  A)4  -  g4)2 

=  (/oA_2  +  /(A_] +/2/Iq  +  /3A] +/4/»2-g2)  (80) 

+  (/2A_2  +  f3h_i  +  f4h0  +  f5hl+  f6h2  -  g4)2 

Again,  the  other  three  terms  which  seem  to  involve /3  (convolutions  centered  at /j,/3,  and/5)  are  not  in  the 
summation.  The  derivative  of  //„  with  respect  to /3  is  derived  in  Eq.  (81),  and  generalized  in  Eq.  (82). 

dH 

*rr  =  2((/®*)2-S2)^1+2((/®/i)4-S4)A-1  (81) 

aJ  3 

^  =  0  .  *2  +  ((/  ®  h)2  -  g2)/t J  +  0  •  h0  +  ((/  ®  A)4  -  g4)*-i  +  0  •  a_2  (g2) 

=  m®h)-g)®hrev)  3 

Again,  we  compute  the  derivative  at  all  points  where  g ,•  exists  (that  is,  at  g2  and  g4). 

We  summarize  the  computation  of  the  derivative  of  noise  term  with  respect  to  a  certain  pixel  in  the  following 
steps: 

1)  From/, compute  fh  =  f®h\ 

2)  Compute  n  =  {ni\ni  =  5(0  •  {fh~g)i\ ,  where  5(0  is  1  only  at  those  points  where  gt  has  a  measurement, 
and  0  otherwise; 


3)  Compute  —  =  (n<8>  hrev). . 

Step  3  can  be  clarified  with  the  illustration  in  Fig.  61,  where  is  non-zero  only  when  i  is  even  in  this  exam- 
ple. 

«0  ”2  «4  "6  "8 


^2  \  h.\ 

\  \  K  h-l  h-2 


kernel  used  to  find- 


kernel  used  to  find 


kernel  used  to  find 


Figure  61.  The  reverse  convolution  in  the  derivation  of  the  noise  term. 
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We  observe  that  the  same  kernel,  hrev,  is  actually  used  every  time,  but  since  some  elements  of  n  (those  with 
odd  indices)  are  zero,  speed  up  may  be  accomplished  by  only  using  the  elements  of  hrev  corresponding  to 
even  indices  of  n. 

The  convolution  of  2-D  images  with  the  reverse  kernel  follows  the  same  rule  as  that  under  1-D  case.  Take  the 
2-D  example  that  only  one  column  of  data  (the  fourth  column)  are  missed,  and  the  size  of  the  blur  kernel  (h) 
is  3  x  3,  as  illustrated  in  Fig.  62.  Similarly, /y  represents  a  pixel  from  the  original  image,  gy  represents  a  pixel 
from  the  measured  image,  the  question  mark  (?)  represents  unmeasured  gy  at  that  pixel,  and  h  is  the  3  x  3 
blur  kernel. 


Figure  62.  A  2-D  example  with  one  column  of  data  (the  fourth  column)  missed. 


Similar  to  the  derivation  steps  in  1-D  case,  we  first  write  out  all  the  terms  involving  a  pixel  (f2 2)  at  which  a 
measurement  (£22) was  made  in  the  noise  term  Hn  as  Eq.  (83), 


e22  =  ((/®/»)1I-g„)2  +  ((/®/»)12-g12)2  +  ((/®/«)21-«21)2 

+  ((/  ®  h)22  -  g22)2  +  ((/  ®  fc)31  -  g31)2  +  ((/  ®  h)32  -  g32)2 


where  the  convolution  of/with  the  blur  kernel  h  is  extended  as  Eq.  (84). 

1  1 

if®h)tj=  X  I  fi  +  m,j  +  n  K,n  (84) 

m  =  -1  n  -  -1 

The  derivative  of  Hn  with  respect  to  pixel /22  can  then  be  derived  as  Eq.  (85),  and  further  generalized  as  Eq. 

(86), 


a-r-  =  2((/<8>/011-g11)-/i1  j  +  2((/<S) /i)12-g12)  •  /ij  0  +  0  ■  /ij 

22 

+  2 ((/  ®  h)2 1  -  g2\ )  •  hQ 1  +  2((/  ®  h)22  -  g22)  •  Hq  q  +  0  ■  /i0j 
+  2((/  <8)  /z)31  —  g3j )  •  h_ j  j  +2 ((/  ®  /i)32  -  £32)  •  h-\to  +  0  '  ^-1,-1 


(85) 


22 


m®h)~8)®hrev)22 


(86) 


where  hrev  has  the  form  as  indicated  in  Eq.  (87),  and  ((/  ®  h)  -  g)  is  computed  at  all  points  where  gy  exists 
(that  is,  at  gu,  g12,  g21,  g22,  gu>  and  g32). 

Ai,i  ^1,0  h\,-\ 
hrev  =  A0  j  ^0, 0  ^0,-1 
^-1,1  0 

We  also  write  out  all  the  terms  involving  an  unmeasured  pixel  (say  /23)  in  the  noise  term  Hn  as  Eq.  (88). 
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(88) 


£23  =  ((/®/»),2-«12)2  +  ((/®/')l4-«14)2  +  ((/®/l)22-S22)2 
+  ((/  ®  h)u  -  g24)2  +  ((/  ®  h)n  -  gn)2  +  ((/  ®  A)^  -  g34)2 

The  derivative  of  Hn  with  respect  to/23  derived  in  Eq.  (89),  and  generalized  in  Eq.  (90). 

*jrr-  =  2((/®«i2-Si2)’*i  j+0-A,  0  +  2((/  ®  «14-*i4)  ■  Alf_, 

0/23 

+  2((/  ®  ^)22  ~  £22) '  V  1  +  ®  *  Vo  +  2((/®  ^)24-g24)  •  A0,-i 

+  2((/  ®  /z)32  -  g32)  •  j  +  0  •  0  +  2((/  ®  A) 34  -  g34)  ■  _i 

(((/®«-«)®*w)2, 

o/23 


(89) 


(90) 


where  hrev  has  the  same  form  as  Eq.  (87),  and  ((/  ®  h)  -  g)  is  computed  at  all  points  where  exists  (that  is, 
at  £l2>  8 \4>  822 *  £24*  £32*  and  £34)* 


We  summarize  the  computation  of  the  derivative  of  noise  term  with  respect  to  a  certain  pixel  under  2-D  case 
in  the  following  steps: 


1)  From/,  compute  fh  —  /  ®  h ; 


2)  Compute  «  =  =  8(i,  7)  •  (/,,  -  g)/;.}  ,  where  8(z,  j)  is  1  only  at  those  points  where  g(J  has  a  measure¬ 

ment,  and  0  otherwise; 


3)  Compute  j- A  =  (n  ®  hrev)ij . 
aJ  ij 

The  following  program  implements  this  algorithm  in  2-D. 

for  (r=0;  r<rows;  r++) 
for  (c=0;  c<cols;  C++) 

fh[r][c]  =  convolve(f,  h,  r,  c); 
for  (r=0;  r<rows;  r++) 
for  (c=0;  c<cols;  C++) 
if  (measured(r,c)) 

n[r][c]  =  fh[r][c]-g[r][c]; 

else 

n[r][c]  =  0; 

for  (r=0;  r<rows;  r++) 
for  (c=0;  cccols;  C++) 

gradf[r]fc]  =  convolve_speedup(n,  hrev,  r,  c); 


7.11  Complexity  Analysis 

Before  comparing  performance  of  the  consistency  methods  and  the  MFA  method  using  experimental  tech¬ 
niques,  we  first  analyze  the  algorithm  complexity.  We  count  each  addition,  subtraction,  multiplication,  divi¬ 
sion,  or  square  root  as  one  flop. 

For  the  consistency  method  using  NLSE  criterion,  if  the  dimension  of  image  ismxn  where  m  and  n  are 
large,  and  b  is  the  number  of  brightness  level,  the  restoration  work  is  dominated  by  the  operations  of  HHQR 
for  the  two  linear  systems  and  the  consistency  testing  which  takes  b  iterations  for  each  row.  The  complexity 

of  HHQR  algorithm  is  on  the  order  of  0^2mn2-  [Trefethen  and  Bau  97],  and  that  of  the  consistency 
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testing  is  on  the  order  of  0(2bmn2) .  Therefore,  the  complexity  of  the  consistency  method  is  on  the  order  of 
0^{2  +  2 b)mn2  -  j ,  where  two  HHQR  is  counted. 

The  complexity  for  the  MFA  method  is  not  deterministic  since  the  times  of  annealing  is  image-dependent.  It 
takes  0((2pq  +  5)mn)  flops  for  noise  term  related  operations  0({2s  +  2t  +  2 st  +  32)mn)  flops  for  prior  term 
related  operations,  and  0(5mn)  for  the  update  operations,  where  p  x  q  is  the  dimension  of  blur  kernel,  s  x  t 
is  the  dimension  of  convolution  kernel  that  simulates  nth  derivatives.  Thus,  the  complexity  of  each  annealing 
is  on  the  order  of  0(( 42  +  2 pq  +  2s  +  2t  +  2 st)mn) .  The  number  of  annealing  iterations  is  determined  by  the 
initial  temperature,  the  final  temperature,  and  the  decreasing  step  of  the  temperature,  which  are  chosen 
image-dependent.  If  the  initial  temperature  is  10,  final  temperature  0.1,  and  decreasing  step  is  0.99,  then  the 
annealing  times  is  10000. 

Therefore,  generally,  for  image  dimension  on  the  order  of  1000  x  1000,  the  MFA  method  has  higher  com¬ 
plexity  than  the  consistency  method.  However,  there  is  a  straightforward  hardware  implementation  of  MFA 
restoration  [Bilbro  and  Snyder  90][Bilbro  et  al.  98]  which  is  very  efficient. 


July  15,  1999 


72 


